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Allocating  and  Launching  MaRV  Interceptors 
Based  on  Miss  Distance  Predictions 


G.  C.  Corynen 

Abstract  and  Summary 

We  have  developed  battle  management  algorithms  for  assigning  and  launching  ground- 
based  interceptors  to  neutralize  a  threat  consisting  of  multiple  decoyed  MaRVs  and  BRVs 
in  their  terminal  phase. 

Decisions  are  based  on  the  minimization  of  expected  leakage  risk,  and  they  specifically 
account  for  uncertainties  and  errors  in  the  identity  of  each  threat  object,  although  detailed 
discrimination  issues  are  not  discussed.  In  contrast  to  standard  probabilistic  approaches, 
the  algorithms  are  physics-based  and  use  kinematic  and  geometric  models  for  estimating 
interceptor-MaRV  miss  distances  and  probabilities  of  kill  to  determine  optimal  waiting 
times  before  launch. 

Based  on  a  test  comparing  interceptor  and  MaRV  technologies,  the  minimax  algorithm 
finds  the  waiting  times  for  which  the  minimum  probability  of  kill  is  a  maximum.  To 
avoid  a  full  multidimensional  optimization,  we  use  a  “greedy”  one-dimensional  approach 
where  only  the  most  urgent  threat  is  addressed  optimally,  conditional  to  allocating  to  the 
remaining  threats  (dimensions)  a  nominal  quantity  of  interceptors  determined  by  previous 
optimizations.  One-on-one  engagements  are  thus  emphasized,  and  the  preponderance  of 
the  work  discussed  in  this  report  was  directed  towards  the  development  of  a  miss  distance 
prediction  algorithm,  called  the  MD  Algorithm. 

While  functional  and  software  testing  of  the  whole  battle  manager  is  still  under  way, 
experiments  with  the  MD  algorithm  have  been  completed  and  have  provided  useful  insights 
which  generally  agree  with  conclusions  reached  by  others.  It  was  found  that  miss  distances 
are  not  very  sensitive  to  relative  interceptor-MaRV  velocities  but  are  quite  sensitive  to 
relative  maneuvering  capabilities,  as  measured  by  relative  lateral  acceleration  limits.  More 
specifically,  for  MaRVs  and  interceptors  with  conventional  capabilities,  and  for  a  total 
defense  system  response  time  of  0.1  second,  the  minimax  miss  distance  was  predicted  to 
be  about  600  feet. 
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Chapter  1 

INTRODUCTION 

Recent  technological  advances  in  the  design  of  Ballistic  Re-Entry  Vehicles  (BRVs)  have 
stimulated  an  increasing  interest  in  the  effectiveness  of  Maneuverable  Re-Entry  Vehicles 
(MaRVs),  and  in  the  development  and  utilization  of  vehicles  designed  to  intercept  MaRVs. 
The  research  documented  in  this  report  is  part  of  a  comprehensive  effort  to  develop  model¬ 
ing  and  simulation  software  to  analyze  terminal  engagements  between  interceptor  defense 
sites  and  multiple  decoyed  MaRVs  and  BRVs.  The  development  of  this  capability  requires 
the  solution  to  various  interceptor  resource  allocation  and  launch  logic  problems,  and  in 
this  report  we  address  two  such  battle  management  problems: 

1.  The  assignment  of  interceptors  to  MaRVs  and  BRVs, 

2.  The  selection  of  interceptor  launch  time  and  direction. 

Significant  contributions  have  been  made  by  others  in  the  development  of  strategies 
to  manage  terminal  engagements  with  MaRVs  (see  References  [1-7],  for  instance).  But 
these  and  other  current  approaches  are  mostly  probabilistic  and  do  not  sufficiently  reflect 
the  physics  of  engagements,  hence  they  are  often  too  abstract.  In  contrast,  the  approach 
presented  here  uses  kinematic  and  geometric  engagement  models  to  predict  interceptor  miss 
distances  for  different  launch  delay  times.  By  varying  these  delay  times,  miss  distances 
can  be  minimized  and  kill  probabilities  maximized.  This  results  in  optimal  interceptor 
resource  allocations,  and  in  the  minimization  of  defense  leakage  risk. 

Although  we  impose  no  specific  restrictions  on  the  quantity  of  interceptors  and  attack¬ 
ing  objects,  we  have  emphasized  one-on-one  engagements,  for  two  major  reasons.  First, 
collective  interactions  between  objects  turn  out  to  be  rather  weak  compared  with  indi¬ 
vidual  interactions  between  one  interceptor  and  one  attacking  object,  particularly  when 
that  object  is  a  MaRV.  Second,  these  one-on-one  interactions  also  turn  out  to  be  relatively 
more  complex  and  computationally  challenging  than  collective  ones,  at  least  in  our  current 
battle  management  framework. 

Similar  to  judging  the  performance  of  a  defensive  system  by  the  proportion  of  threats 
which  leak  through  it,  the  performance  of  an  interceptor,  and  its  associated  launch  policy, 
may  be  judged  by  the  probability  that  the  interceptor  will  neutralize  the  threat  to  which  it 
was  assigned.  Because  engagements  typically  take  place  a  considerable  time  after  launch, 
the  battle  manager  must  predict  this  probability,  sometimes  far  in  advance,  to  determine 
when  to  launch,  and  to  decide  how  many  interceptors  should  be  assigned  to  the  threat. 
Considering  that  the  kill  probability  of  an  interceptor  is  very  sensitive  to  the  miss  distance 
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achieved  by  it,  a  capability  to  predict  miss  distances  is  essential  to  an  effective  battle 
manager.  This  prediction  problem  is  the  principal  subject  addressed  in  this  report. 

Our  discussion  is  structured  as  follows.  In  the  remainder  of  this  chapter,  we  summarize 
the  battle  managment  problem  and  present  our  general  solution  to  it.  In  particular,  we 
solve  the  first  problem  listed  earlier:  that  of  assigning  interceptors  to  MaRVs  and  BRVs. 
Subsequent  chapters  are  devoted  to  solving  the  second  problem  of  selecting  the  time  and 
direction  of  interceptor  launch.  In  Chapter  2,  we  present  the  mathematical  preliminar¬ 
ies  required  to  discuss  this  second  problem.  In  Chapter  3,  we  present  our  miss  distance 
prediction  and  minimization  methodology,  and  in  Chapter  4  we  report  the  result  of  com¬ 
putational  experiments  with  the  miss  distance  algorithm.  Conclusions  and  future  work 
recommendations  are  made  in  Chapter  5.  Nine  comprehensive  appendices  round  out  the 
report. 


1.1.  The  Battle  Management  Problem 

A  simplified  pictorial  of  the  general  battle  management  problem  considered  in  this 
paper  is  described  in  Figure  1.1,  where  a  set  of  threat  objects  O  =  {Og:  q  —  1,  . . . ,  Q}  is 
threatening  a  collection  of  assets  A  =  {a*.:  k  =  1,  . . . ,  K}  which  are  protected  by  pools  of 
interceptors  J  =  { Iu :  u  =  1,  . . . ,  U}  located  at  different  sites.  Threat  objects  are  of  two 
basic  types:  MaRVs  and  BRVs.  Interceptors  are  assets  which  can  be  members  of  several 
pools.  They  are  also  of  two  basic  types:  Overlay  Interceptors  (type  1),  and  Underlay 
Interceptors  (type  2). 

As  indicated  earlier,  the  major  objective  of  the  battle  manager  is  to  make  the  best  use 
of  available  resources  to  defend  its  assigned  assets.  This  is  a  problem  in  optimal  decision 
making  which  involves  three  principal  subproblems,  as  shown  in  Figure  1.2.  In  practice, 
the  battle  manager  rarely  knows  with  certainty  that  the  incoming  objects  or  threats  are 
in  fact  MaRVs  and,  in  the  next  sections,  we  discuss  how  these  types  of  uncertainty  are 
incorporated  into  our  framework.  But  we  do  not  discuss  detailed  discrimination  issues. 

Referring  to  Figure  1.2,  threats  must  first  be  prioritized  so  that  “important”  threats 
can  be  processed  first.  This  prioritization  is  particularly  important  in  threat-rich  situations 
where  only  a  subset  of  the  threats  can  be  processed  by  the  battle  management  computers 
during  one  update  cycle,  and  we  discuss  it  in  greater  detail  in  Section  1.3. 

Next,  the  battle  manager  must  decide  how  many  interceptors  should  be  allocated  to 
each  incoming  threat,  and  which  type  of  interceptor  should  be  used.  In  some  schemes,  this 
is  accomplished  simply  by  specifying  a  desired  kill  probability,  and  allocating  a  sufficient 
number  of  interceptors  to  achieve  this  specification.  In  our  scheme,  we  also  allow  this 
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Figure  1.1.  Typical  engagement  geometry. 


approach,  but  generally  take  the  broader  view  that  such  decisions  involve  several  additional 
parameters  relating  to  defense  attrition,  regret,  and  so  on.  This  is  discussed  further  in 
Section  1.2  below. 
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Figure  1.2.  The  battle  manager  performs  three  major  functions. 


In  the  interceptor-launch  logic  subsystem,  the  battle  manager  decides  the  direction 
and  time  of  interceptor  launch.  This  is  the  problem  emphasized  in  this  report,  and  we 
discuss  it  in  Chapters  2  through  4. 


1.2.  Battle  Management  Decision  Policy 

In  assigning  interceptors  to  threat  objects,  the  battle  manager  must  achieve  an  ap¬ 
propriate  balance  between  the  costs  resulting  from  threat  leakage  and  the  expenditures 
required  to  prevent  such  leakage.  Many  of  the  parameters  determining  these  costs  and 
expenditures  are  imperfectly  known.  Leakage  costs,  for  instance,  depend  on  the  type  and 
size  of  the  threat  and  on  the  vulnerability  and  lethality  of  threat  objects  and  defended 
assets.  Defense  costs  are  dominated  by  the  intrinsic  (“lifecycle”)  costs  of  interceptors  and 
by  the  risk  of  reduced  protection  associated  with  the  launch  of  interceptors  from  a  limited 
inventory. 

We  approach  the  battle  management  problem  as  a  stochastic  optimization  problem 
where  an  optimal  balance  is  sought  between  leakage  risk  and  defense  cost.  For  tills  opti¬ 
mization,  we  use  total  expected  loss  L  as  the  objective  function.  More  specifically,  our  goal 
is  to  minimize  the  loss  function 

Q 

Lr(maq)  =  ^2(Lq'PLl(maq)  +  Cq(maq))  ,  (1.1) 

q=l 
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is  the  loss  incurred  if  object  Oq  leaks  through, 


where  Lg 

pLq  is  the  probability  that  0q  leaks  through, 

Cq  are  the  costs  incurred  by  allocating  defense  resources  to  0q, 
maq  is  the  quantity  of  defense  resources  of  type  a  allocated  to  0q. 

This  is  a  Q-dimensional  constrained  optimization  problem  where  an  optimal  vector 
771*  =  [m*j ,  m*9,  777*^]  must  be  found,  subject  to  the  constraint  that  defense 

resources  of  type  a  cannot  exceed  Ua : 

Q 

^2maq<Ua  .  (1.2) 

9= 1 

Although  this  chapter  is  only  a  summary  of  the  battle  management  approach,  a  few  com¬ 
ments  should  be  made  to  demonstrate  the  generality  and  flexibility  of  the  objective  function 
in  Equation  (1.1). 

Even  though  the  additive  structure  of  this  function  might  suggest  that  the  probabilistic 
events  and  random  variables  involved  should  be  either  disjoint  (mutually  exclusive)  or 
independent,  this  is  not  the  case,  as  shown  in  Reference  [8],  p.  12.  Additivity  or  linear 
“accumulation”  of  the  losses  Lq  are  also  not  required.  All  we  need  is  to  assign  or  to  estimate 
a  loss  Lq  for  each  threat  object  Oq.  Admittedly,  if  several  0qs  are  headed  for  the  same  asset 
(target),  these  loss  values  may  depend  upon  which  0q  reaches  its  destination  first,  but  that 
is  an  assessment  problem,  not  a  problem  in  structuring  the  loss  function  of  Equation  (1.1). 
Stepping  away  from  the  mathematics  for  a  moment,  the  battle  manager  assesses  these 
Lq  values  in  real  time,  based  on  radar  and  other  sensor  updates  concerning  the  position, 
direction,  and  destination  of  the  associated  0q s.  If  several  objects  are  perceived  by  the 
sensors  to  be  headed  for  the  same  target,  this  information  will  be  used  to  estimate  which 
object  will  reach  the  target  first,  the  damage  caused  by  that  object,  and  any  additional 
damage  which  will  be  caused  by  the  remaining  objects  still  in  flight. 

Whereas  maq  can  be  used  to  indicate  many  different  types  of  resources  such  as  warhead 
yield,  booster  capability,  payload,  and  so  on,  we  shall  use  it  exclusively  to  specify  the 
quantity  of  interceptors  of  type  a  allocated  to  threat  0q. 

Turning  our  attention  to  the  computation  of  Equation  (1.1),  we  shall  see  that  pLq  and 
Cq  are  nonlinear  functions  of  maq,  and  these  parameters  also  share  common  variables  with 
Lq.  For  instance,  all  three  parameters  depend  on  whether  threat  object  0q  is  a  MaRV  or  a 
BRV.  Hence  the  optimization  problem  cannot  easily  be  twisted  into  a  linear  programming 
problem,  unfortunately.  Considering  the  computational  implications  of  attempting  to  solve 
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the  full  nonlinear  problem,  we  have  chosen  a  more  myopic  or  “greedy”  approach.  Rather 
than  solving  the  full  Q-dimensional  problem,  we  consider  one  threat  object  Oq  at  a  time, 
starting  with  the  highest  priority  one,  and  we  attempt  to  find  the  optimal  value  m*g  for 
that  threat,  invoking  some  mild  assumptions  about  the  future  treatment  of  the  remaining 
Q  —  1  threat  objects.  More  specifically,  we  shall  assume  that  BRVs  down  the  priority  list 
are  assigned  a  single  interceptor  of  the  appropriate  type,  and  that  MaRV’s  are  assigned 
two,  one  of  each  type,  until  the  limit  Ua  is  reached.  This  assumption  sets  a  “level”  for  the 
local  optimization  where  the  optimal  value  m*aq  for  the  highest-priority  Oq  is  sought.  This 
one-dimensional  approach  obviously  encourages  iteration  over  different  levels.  If  we  find, 
for  instance,  that  m*9  consistently  exceeds  2,  the  nominal  allocation  level  of  the  remaining 
objects  should  probably  be  increased,  but  this  can  easily  be  determined  experimentally  on 
the  computer.  In  the  worst  case,  this  “local”  scheme  may  not  approach  the  true  optimum— 
as  obtained  by  Q-dimensional  optimization — and  a  more  sophisticated  approach  may  be 
necessary.  This  will  be  established  with  exhaustive  tests  when  the  battle  manager  is  fully 
implemented. 

In  the  next  two  sections,  we  discuss  each  of  two  components  of  the  objective  function 
in  greater  detail. 

1.2.1.  COMPUTING  LEAKAGE  LOSSES 

We  estimate  the  costs  of  threat  leakage  by  the  expected  asset  loss  resulting  from  the 
leakage.  If  Lq  is  the  loss  incurred,  conditional  to  object  Oq  leaking  through,  and  pLq  is  the 
probability  that  Oq  leaks  through,  the  expected  leakage  loss  is 

<? 

L  =  L1  ’PLq(mocq)  ,  (!-3) 

where  we  emphasize  the  fact  that  the  pLq  depends  on  the  quantity  maq  of  interceptors  of 
type  a  assigned  to  threat  0q  by  the  battle  manager.  While  L  may  include  interceptor  losses 
due  to  leakage,  it  does  not  account  for  the  loss  of  protection  arising  from  the  assignment 
and  launch  of  interceptors.  Such  losses  will  be  discussed  below,  and  in  Section  1.2.2. 

But  Lq  and  pLq  both  depend  on  the  type  of  the  threat  object  Oq.  This  introduces 
dependencies  between  pLq  and  Lq,  and  requires  the  modification 

Q 

L  =  [(L<lPLq\oq  =  M)  Pq1  +  (LqPLq\oq=B)  Pq] 

q=l 


q=l 
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where 


+Pf+Pq  =  h 

is  the  probability  that  0q  is  a  MaRV, 
pq  is  probability  that  0q  is  a  BRV, 
pjj*  is  the  probability  that  0q  is  a  decoy, 

L^(Lg)  is  the  leakage  loss  if  0q  leaks  through  and  is  a  MaRV  (BRV), 
pM{P»q)  is  the  leakage  probability  if  0q  is  a  MaRV  (BRV). 


But  L M  and  depend  on  the  probable  destination  and  reach  of  0q,  and  on  the  assets 
threatened.  Hence 

K 

=  E  vMp"  (1-5) 

k=  1 

and 

K 

if  =  E  .  (1-6) 

k=l 

where  V(aj.)  is  the  value  of  asset  a*., 

pM  is  the  probability  that  0q  is  headed  for  a j.  if  0q  is  a  MaRV, 
and  similarly  if  0q  is  a  BRV. 

In  our  framework,  we  estimate  or  approximate  p(jk  with  the  probability  that  asset  lies  in 
the  footprint  fq  of  0q  (see  Reference  [9]  for  a  discussion  of  footprint  calculations).  Thus, 

Pql  =  Pq  €  fq  Io9=m)  (1-7) 

and 

Pqk  =  Pq  ( ak  £  fq\oq  =  B )  •  (1*®) 

where  fq  is  the  footprint  of  0q  (a  random  function). 


Next,  considering  maq  interceptors  assigned  to  a  single  0q,  and  assuming  that  the 
interceptors  damage  0q  independently, 


and 

I 


(1.9) 

(1.10) 
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Using  the  simple  facts  that 

PLq(m<*q  =  1)  Q^=B  =  (i  “  =  !)  (L11) 

and 

= i)|0]=M = u  -  **.<">«. = d|0i=m)  -  <L12> 

where  pKg(maq  =  1)  is  the  probability  that  an  interceptor  of  type  a  will  kill  0q,  we  can 
combine  Equations  (1.1) — (1.12)  to  obtain  the  overall  leakage  loss 


where 


and 


=  (*  -  Pk [  (m*9  =  !))  a*  » 

Q  is  the  quantity  of  threat  objects, 
K  is  the  quantity  of  assets, 


(1.13) 

(1.14) 


maq  is  the  quantity  of  interceptors  of  type  a  launched  simultaneously 
towards  0q, 

Pjcq(motq  =  1)  is  the  single-shot  kill  probability  for  a  MaRV  threat  0q, 
and  similarly  for  the  BRV  superscript  B. 


The  individual  kill  probabilities  p and  p^  used  in  Equations  (1.13)  and  (1.14)  may 
require  minor  corrections  in  cases  where  several  interceptors  of  a  given  type  are  launched 
from  different  sites  towards  a  given  threat  object.  These  probabilities  do  in  fact  depend 
upon  the  site  location  from  which  the  interceptors  are  launched.  Sometimes,  a  preferred 
site  may  be  depleted,  and  may  request  the  assistance  of  other  sites  to  achieve  the  total 
number  maq  of  interceptors  required.  If  maqs  is  the  number  launched  from  site  s, 

5 

ocq  —  ^  ^  'm&qs  i  (1*1^) 

where  S  is  the  total  number  of  sites. 


Referring  to  Equation  (1.9),  another  expression  must  be  substituted  for  the  probability 
pL  ( • )  in  accordance  with  the  equality 


\Oa  =  M 


ft  \  PLgs(mc 


(1.16) 
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where  PLqs{  * )  is  the  probability  that  0q  will  leak  through  if  a  single  interceptor  is 
launched  from  site  s. 

But  finding  the  optimal  values  for  maqs  itself  requires  a  multidimensional  optimization.  To 
avoid  such  computations  in  the  calculation  of  the  optimal  number  m*aq ,  we  assume  that  in¬ 
terceptors  are  launched  from  the  most  effective  site.  In  these  calculations,  Equation  (1.16) 
is  replaced  by  the  equation 

Wag 

PLg(m«g)b9=M  =  II  “gj,  {PLqs(mogs  =  1)Io,=m}  >  (L17) 

u—  1 

where  Sau  is  the  set  of  sites  at  which  at  least  one  interceptor  of  type  a  is  available  after 
the  nth  interceptor  has  been  allocated  (similarly  when  Oq  =  B). 

If,  instead  of  the  two-tiered  approach  where  only  interceptors  of  the  same  type  are 
allowed  in  any  given  allocation,  we  allow  a  mix  of  types,  then  the  above  minimization  can 
also  be  extended  over  types.  In  such  cases,  we  obtain  a  probability  which  is  independent 
of  type: 

mq 

PLq(mq)\oq=M  =  II  =1)lo9=M/  •  (L18) 

u=  1  aE{  1,2} 

The  effects  of  depleting  interceptor  resources  can  be  studied  very  easily  with  this  model. 
Whenever  a  site  runs  out  of  interceptors,  the  total  set  of  nonempty  sites  is  reduced,  and 
the  minimization  is  bound  to  yield  a  larger  number.  In  the  extreme,  there  may  not  be 
any  interceptors  left  at  all,  and  pLg(mq)  must  be  set  to  unity.  The  effects  of  varying 
the  resources  allocated  to  each  site  can  thus  be  studied  in  order  to  find  an  appropriate 
distribution  of  resources  over  defense  sites. 

For  such  studies,  the  loss  function,  Equation  (1.13),  may  be  decomposed  into  two  parts, 
one  indicating  no  depletion  and  the  other  total  depletion,  as  follows: 

l  =  £  £  h°*>  ) 

q~l  k= 1 

+  £  ^2[V(ak)(p^kPf  +PqkPq)]  >  (1’19) 

q=Qi+ 1  fc=i 

where  Oqt  is  the  last  threat  object  on  the  priority  list  for  which  a  sufficient  quantity  of 
interceptors  is  available. 


10 


1.2.2.  COMPUTING  DEFENSE  COSTS 


In  our  simplified  models,  the  defense  cost  variable  Cq  of  Equation  (1.1)  represents  only 
the  intrinsic  (“lifecycle”)  costs  of  the  interceptors,  and  has  the  simple  form 

Cq  =  maqC(Ia)  ,  (1.20) 

where  C(Ia)  is  the  lifecycle  cost  of  one  interceptor  of  type  a. 

An  earlier  version  of  our  model  included  in  Cq  various  “loss  of  protection”  costs  result¬ 
ing  from  the  depletion  of  interceptors.  These  costs  have  now  been  included  in  the  overall 
risk  function  of  Equation  (1.19). 

As  a  concluding  comment,  Equation  (1.13)  is  a  convenient  form  for  L  because  we 
can  minimize  L  by  maximizing  each  pKq,  q  =  1,  ...,  Q  (other  things  held  fixed).  The 
performance  of  the  defense  in  engaging  Q  threats  can  then  be  judged  by  how  well  each 
individual  threat  is  engaged,  and  that  is  why  we  have  emphasized  the  one-on-one  situation 
in  the  remainder  of  this  report. 

1.3.  Threat  Prioritization 

To  avoid  the  saturation  of  its  sensors,  radars,  and  computers  when  the  defense  is  con¬ 
fronting  massive  MaRV  attacks,  the  battle  manager  processes  and  updates  its  information 
only  for  a  critical  subset  of  threat  objects.  This  subset  is  determined  by  an  appropriate 
balance  between  the  update  rate  (clock  time  increment)  and  the  importance  of  threats. 
Three  criteria  are  combined  into  a  threat  ranking  rule: 

1 .  severity, 

2.  vulnerability, 

3.  urgency. 

1.3.1.  THREAT  SEVERITY 

Used  to  reflect  the  damage  expected  from  the  leakage  of  a  threat  object  if  it  leaks 
through,  this  criterion  is  simply  the  conditional  expected  asset  loss  Lq  derived  earlier 
(Equations  (1.5),  (1.6),  (1.13)): 

k= 1 
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1.3.2.  THREAT  VULNERABILITY 


There  is  no  point  in  launching  an  interceptor  if  its  pK  is  zero,  and  we  include  in  our 
ranking  rule  the  probability  p°Kq  that  0q  will  be  killed  if  an  interceptor  is  launched  without 
delay,  be.  when  twait  =  0.  Target  Oq  is  thus  ranked  lower  than  target  Or  if 

it  =  0)  <  Ptfrityiait  —  0)  •  (1.22) 

As  usual,  p°K  depends  on  the  type  of  Oq ,  and  an  expected  value  is  used: 

rl,  =  (?*, Io,=m)  p,"  +  (p°Jo,=b)  •  <L23> 

The  remaining  chapters  in  this  report  describe  how  p®  is  predicted. 

1.3.3.  THREAT  URGENCY 


When  an  interceptor  launch  is  delayed  too  long,  the  interceptor  (I)  may  fail  to  intercept 
its  assigned  threat  object  (Oq),  or  it  may  intercept  it  inside  the  keep-out  zone,  which  is 
not  allowed.  The  maximum  amount  of  waiting  time  allowed  is  thus  a  good  measure  of  the 
urgency  of  a  threat,  and  we  use  that  measure  in  our  ranking  criterion.  This  waiting  time 
(rg)  is  simply  the  difference  between  the  time  required  by  the  threat  object  to  reach  the 
keep-out  zone  and  the  time  required  by  the  interceptor  to  reach  the  keep-out  zone,  at  the 
point  Xko •  With  minor  modifications  due  to  the  differences  between  the  predicted  impact 
point  and  the  predicted  keep-out  pierce  point,  we  can  use  the  results  of  Section  2. 2. 1.1  to 
obtain  rq. 


If  (t[o)  are  the  times  required  by  Oq  (I)  to  reach  Yfc0,  then 


_  _  ±M 
r?  —  *ko 


lko 


(1.24) 


Using  Equation  (2.58)  to  solve  for  the  interceptor  flight  time  t|o,  and  Equations  (2.10)- 
(2.13)  with  the  required  modifications  to  account  for  the  keep-out  zone,  we  obtain 


where 


iko  =  max{0,  r0}  +  min{^',  tF} 


F  4F  i 


Tn  = 


Hq  —  X®  •  ey 

eV$'ey 


iF 

*0 


_  -2(gp  -  hko ) 

^(in°l  +  inFIKo-e,  ’ 


t *  = 


-2  (X°-ey-hko) 
(\Vq°\  +  WF\)e^ey  ’ 


(1.25) 

(1.26) 

(1.27) 

(1.28) 


12 


hka  is  the  height  of  the  keep-out  zone, 

Ho  is  the  nominal  height  of  the  atmosphere, 

X®  is  the  current  position  of  Oq, 

Vq  is  the  current  velocity  of  Oq , 
ey  is  a  unit  vector  along  the  altitude  coordinate  y, 
eyo  is  a  unit  vector  along  the  current  velocity  vector  V9°, 
Vq  is  the  terminal  velocity  of  0q. 

These  results  will  be  derived  and  explained  in  Section  2.2. 1.1. 
1.3.4.  OVERALL  THREAT  RANKING  CRITERION 


To  obtain  the  overall  ranking  rule,  the  three  criteria  derived  above  are  heuristically 
combined  into  the  simple  ratio 


(1.29) 


Observe  that  the  optimal  probability  p*Kq  (and  also  depends  on  the  location  of  the 
interceptor.  Since  interceptors  can  be  launched  from  several  sites,  we  account  for  this 
dependence  by  defining  p*  as  the  maximum  probability  over  all  potential  launch  sites,  as 
derived  in  Equation  (1.17). 


1.4.  Interceptor  Allocation 

The  principal  task  in  allocating  interceptors  is  to  determine  the  optimal  number  m*aq 
of  type  a  interceptors  to  be  assigned  to  the  highest-priority  threat  object  Oq.  As  discussed 
earlier,  we  approach  this  problem  as  a  two-tier  optimization  problem  where  only  two 
types  of  interceptors  are  used:  an  overlay  interceptor  (a  =  1)  and  an  underlay  interceptor 
(a  =  2).  Overlays  are  typically  used  at  altitudes  exceeding  100,000  feet,  where  passive 
maneuvering  is  seriously  limited,  and  where  a  thrusting  interceptor  may  have  a  significant 
advantage.  Underlays  are  normally  used  in  a  “last  ditch”  defense  at  lower  altitudes  where 
maneuvering  capabilities  are  essential.  Such  interceptors  typically  have  a  shorter  range, 
but  a  higher  acceleration  capability,  and  our  example  interceptors  in  Chapter  2  reflect  such 
differences.  We  do  not  allow  mixed  launches  in  our  current  version  of  the  battle  manager: 
m*  always  equals  the  total  number  of  interceptors  allocated  to  Oq  in  any  given  allocation. 
Threats  may  reappear  on  the  priority  list  if  they  were  not  killed  by  the  m*aq  interceptors 
allocated  in  a  previous  cycle,  in  which  case  a  new  allocation  is  made,  typically  with  Type  2 
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interceptors.  To  find  ra*?,  we  rise  an  objective  function  which  combines  Equations  (1.13), 
(1.17),  (1.19),  and  (1.20): 


Q 

*  =  £ 
q=l 


fi 

Y  [V (  ak )  )pf  +  PqkP^q  (™«<?  )Pq  )  +  ™*q  <?(  4  ) 

k= 1 


where 


mQq 

=  II  ™ln  =  !)  \oq=M }  > 

c  &au  ^  ' 

11=1 

the  minimum  leakage  probability  achievable  over  all  nonempty  sites. 


(1.30) 


As  indicated  at  the  outset,  we  avoid  the  Q-dimensional  optimization  by  minimizing  L 
over  maq  only  for  q  =  1,  be.,  only  for  the  highest-priority  object  Oq.  The  other  values 
mag ,  q  =  2,  . . . ,  Q,  are  set  to  2  for  MaRVs  and  to  one  for  BRVs.  This  is  only  a  test  level 
for  an  initial  implementation,  and  it  can  obviously  be  refined  if  experiments  indicate  it 
should  be  done.  In  the  worst  case,  no  such  rule  can  be  used,  and  we  are  forced  to  return 
to  the  multidimensional  approach. 

Referring  to  the  sketch  of  Figure  1.3,  the  current  approach  is  to  find  the  minimum 
of  the  one-dimensional  function  Fq  =  Eq  +  Cq  +  W ,  where  Eq  is  the  expectation  of  the 
conditional  loss  Lq  of  Equation  (1.21),  and  q  —  1  since  we  are  considering  only  the  highest 
priority  target.  This  function  is  simply  Equation  (1.30)  with  all  maq  values  fixed,  except 
for  ma i,  and  can  be  expressed  as 

K 

Fi(m.al)  =  Y  V(ak)  (PnP^im^)P^  +  PuPfi(m<>i)P?)  +  m«i °(I«)  +  W  i  C1-31) 

fc=j 

where  W  is  the  (fixed)  risk  associated  with  the  remaining  objects  O2,  . . . ,  Oq.  There  is 
no  closed  form  solution  for  the  minimum  of  this  function,  and  numerical  methods  must 
be  employed.  We  use  Brent’s  algorithm  (Reference  [10])  in  our  Appendix  E  to  find  this 
minimum.  (This  algorithm  will  also  be  used  in  Chapter  3  to  compute  interceptor-MaRV 
miss  distance). 

With  one  simplification,  however,  a  closed  form  solution  to  this  minimization  can 
be  found,  and  we  derived  such  a  solution  for  exploratory  purposes.  If  we  assume,  for 
instance,  that  pLl(ma  1)  is  site-independent,  then  p^(ma  1)  =  (pL1(ma  1  =  l)|0l=^)mal, 
and  similarly  for  0\  —  B. 

If  we  let  w  be  the  number  of  interceptors  left  when  ma  1  =  0,  Equation  (1.31)  can  then 
be  written  as 
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Figure  1.3.  The  optimal  quantity  m*aq  of  interceptors  allocated  to  the  highest-priority 
threat  0q  is  found  by  balancing  leakage  costs  against  defense  costs  sub¬ 
ject  to  the  constraint  that  only  w  interceptors  are  available. 


K 


k= 1 


F\{ma i)  =  y(afc)  (,Pli(™«i  =  1)Io1=m) 


\  min{ty,mai} 


+  PfkPf  (Pii(TOal  =  1)Io1=b) 


+  min{«;,  mQi  }C(Ia )  +  W  .  (1.32) 

To  find  m*  j,  we  find  the  minimum  ra*®  of  Equation  (1.32),  neglecting  interceptor  depletion 
constraints.  Then  we  introduce  this  constraint  by  setting  m*2  =  min{u>,  rn^j}.  To  this 
end,  F\  ( ma i )  has  the  form 


F,(m„,)  =  A(7S)"«+B(if1)»-+ro„,C(/«)  +  W  ,  (1.33) 

^^  =  ^ln(7")(7Sr->>  +  Bln(7®)(7fir«>+C(/a)  .  (1.34) 
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If  Oq  is  within  radar  range,  7^  %  0,  and  if  mai  >  1, 

^fi-AN7SKiar-+c(/.) 

and  the  solution  to  —7^ (mQfl )  =  0  is 

d(mai) 


(1.35) 


711 


*0 

al 


(1.36) 


where  A  =  ^  (ak)pff  •  Reintroducing  the  interceptor  attrition  constraint,  the 

optimal  number  of  allocated  interceptors  is  then 


to”,  =  min{m*0,,  ui} 


(1.37) 


where  w  is  the  number  of  available  a-type  interceptors. 


1.5.  Interceptor  Launch  Logic 

Given  that  the  statistical  outcome  of  an  engagement  is  virtually  determined  by  the 
probability  of  kill  (pK)  of  the  interceptors,  and  that  this  pK  strongly  depends  on  when 
and  where  the  engagement  takes  place,  launch  time  and  launch  direction  are  dominant 
parameters  in  the  management  of  a  battle.  Premature  interceptor  launches  may  cause 
early  burn-out,  giving  a  MaRV  considerable  speed  advantage  at  the  time  of  intercept.  A 
late  launch  may  not  allow  the  interceptor  to  reach  its  maximal  velocity  prior  to  or  during 
the  engagement,  again  limiting  its  effectiveness.  Similarly,  if  an  interceptor  is  launched  in 
an  unfavorable  direction,  a  significant  energy  expenditure  may  be  required  to  force  it  on 
course. 

Consequently,  and  as  shown  in  Figure  1.4,  the  battle  management  launch  logic  subsys¬ 
tem  operates  on  MaRV  state  information  to  produce  two  major  outputs: 

1.  Interceptor  launch  time, 

2.  Interceptor  launch  direction. 

These  decisions  are  made  and  updated  at  a  rate  determined  by  the  rate  at  which  the 
MaRV  state  updates  (Xm)  are  received  from  various  sensors,  and  by  the  internal  clock 
time  of  the  battle  manager.  Depending  upon  the  quantity  of  objects  on  the  threat  priority 
list,  computer  saturation  may  occur,  in  which  case  the  decision  process  must  be  slowed 
down  accordingly. 
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Figure  1.4.  At  each  state  update  time  t,  the  battle  management  launch  logic  sub¬ 
system  operates  on  the  MaRV  state  vector  estimate  Xm^)  to  produce 
the  optimal  launch  direction  A *Xj(t)  and  launch  time  t*L(t). 

The  output  decisions  of  this  logic  subsystem  also  depend  upon  the  technological  ca¬ 
pabilities  and  limitations  attributed  to  both  the  interceptor  and  the  MaRV.  Maneuvering, 
speed,  and  acceleration  capabilities  of  both  vehicles  are  constrained  by  limitations  in  the 
development  and  fabrication  of  metals,  composites,  and  other  materials.  We  represent 
such  constraints  in  our  models  by  defining  limits  on  weight,  lateral  acceleration,  speed, 
and  other  factors,  and  we  represent  these  as  an  input  vector  0  as  shown  in  Figure  1.4. 

The  optimal  launch  time  (and  the  optimal  waiting  time)  is  found  by  maximizing  pK 
over  waiting  time.  Again  pK  is  the  principal  variable,  with  interceptor  lethality  and  miss 
distance  as  major  parameters: 

PKq(rnaq  =  1)  =  prob  {KRaq  >  MDaq }  ,  (1.38) 

where  KRaq  is  the  kill  radius  of  an  interceptor  of  type  a  against  threat  0q, 

MDaq  is  the  miss  distance  achieved  by  an  interceptor  of  type  a  against  0q. 

In  this  report,  we  do  not  further  discuss  lethality  issues,  and  we  concentrate  on  the  miss 
distance  MDaq.  The  sole  purpose  of  Chapter  3  is  to  develop  an  algorithm  for  predicting 
MDaq,  and  for  optimizing  MDaq  and  pKaq  over  waiting  time. 

By  extending  the  MaRV  trajectory  at  the  optimal  launch  time,  we  predict  a  point  of 
intercept,  and  we  use  this  point  to  determine  an  interceptor  launch  direction  vector  A Xj. 

1.6.  Battle  Management  Flowchart 

A  simplified  flowchart  of  the  battle  manager  is  shown  in  Figure  1.5.  First,  threats  are 
prioritized  in  accordance  with  the  criterion  discussed  in  Section  1.3.,  where  kill  probabilities 
are  maximized  over  sites  and  type,  as  derived  in  Equation  1.18  (Block  1).  Then,  the 
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highest-priority  threat  0q  is  placed  on  another  queue  in  accordance  with  its  t*wq  (3).  This 
second  ordering  queue  is  needed  because  highest-priority  threats  are  not  necessarily  the 
most  urgent  ones,  as  measured  by  The  overall  priority  of  an  object  in  this  t*,9-queue 
may  vary  considerably  while  it  is  waiting  to  be  addressed,  and  only  when  t**q  approaches 
the  update  time  increment  r  can  a  final  interceptor  assignment  and  launch  time  be  derived. 
If  there  is  a  threat  on  this  Waiting  queue  whose  /.*,  is  smaller  than  the  update  time,  continue; 
otherwise,  return  to  consider  next  threat  (4). 


Updates:  threat  object  states 


Comments 


Figure  1.5.  Battle  management  flowchart  (continued  on  the  following  page). 
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From  6 


If  yes,  threat  is 
killed.  If  no,  threat 
leaks  through  this 
engagement. 


In  accordance  with 
AiJ,  the  time-to- 
engagement 


If  the  engagement  will  take  place  inside  the  keep-out  zone  for  the  current  threat  on  the 
waiting  list,  then  the  threat  leaks  through,  and  the  assets  within  its  reach  are  destroyed  (5 
and  6).  If  not,  assign  the  optimal  number  of  interceptors  m*  q,  launch  these  inteceptors,  and 
remove  them  from  the  inventory  (8).  Simulate  the  engagement,  obtain  the  simulated  miss 
distance,  and  calculate  the  simulated  (“actual”)  probability  of  kill  (9  and  10).  Using 
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pA  as  the  “true”  probability,  generate  a  number  p  from  a  uniform  random  distribution 
(11),  and  ascertain  whether  to  remove  it  from  the  threat  list.  (14).  If  not,  the  threat  leaks 
through  the  current  engagement,  and  must  be  returned  to  the  active  threat  list  (13),  where 
it  is  processed  as  a  new  threat. 

Future  chapters  will  show  how  some  of  these  computations  are  done.  Chapter  3  in 
particular  will  concentrate  on  the  computation  of  miss  distances  and  waiting  times. 
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Chapter  2 

MATHEMATICAL  PRELIMINARIES 

In  this  chapter  we  present  the  simplified  equations  describing  the  motion  of  MaRVs 
and  interceptors,  and  we  derive  reachability  and  interception  conditions  needed  for  the 
calculation  of  miss  distances.  Both  vehicles  are  treated  as  point  masses,  and  our  principal 
goal  is  to  describe  the  behavior  of  these  vehicles  when  they  are  negotiating  simple  but 
stressing  maneuvers.  First,  we  discuss  general  equations  applying  to  both  types  of  vehicles. 
Then  we  concentrate  on  MaRVs  and  interceptors  individually,  and  we  conclude  with  the 
derivation  of  reachability  and  interception  conditions. 

2.1.  General  Equations 

A  particle  of  mass  m  moving  with  constant  speed  s  and  attracted  toward  a  fixed  point 
O  by  a  force  F  will  travel  along  a  circle  whose  center  is  at  O  and  whose  radius  is  r  =  ms2 /F. 

Consider  next  a  vehicle  such  as  a  MaRV  or  an  interceptor  experiencing  a  constant  lift 
force  (see  Figure  2.1). 


D 


Figure  2.1.  Describing  the  motion  of  a  vehicle  under  constant  lift  conditions. 
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(2.1) 


L  =  \cL(a,M)p\V\2S  , 

where: 

Cl  is  the  lift  coefficient, 
a  is  the  angle  of  attack, 

M  is  the  Mach  number, 

p  is  the  atmospheric  density,  p  =  poe~hlh°, 
ho  —  20,000  feet  (Reference  [11],  p.  F166), 
S  is  the  effective  area, 

|V|  is  the  vehicle  speed. 


For  such  a  vehicle 

m|V|2  _  2m 
r  ~  L  ~  CL( a,  M)pS 


(2.2) 


Observe  that  the  speed  |V|  does  not  appear  in  this  expression,  so  that  we  may  ignore  the 
drag 

D  = -±CDSp\V\2  ,  (2.3) 

where  Cd  is  the  drag  coefficient.  Admittedly,  under  some  conditions,  Cl  and  M  (and 
occasionally  S)  are  themselves  functions  of  velocity,  but  such  conditions  may  be  neglected 
for  the  problem  solved  in  this  report.  For  instance,  the  effects  of  Mach  number  M  on  the 
lift  coefficient  Cl  are  most  pronounced  for  high  values  of  M  (typically  M  >  10),  and  we 
will  show  that  optimal  evasive  engagements  rarely  take  place  in  such  regimes. 


There  is  also  the  possibility  of  “induced  drag,”  where  a  lifting  maneuver  also  causes 
additional  drag,  particularly  for  large  values  of  a.  While  we  are  especially  interested  in 
maximal  maneuvers  where  a  is  in  fact  large,  induced  drag  effects  have  been  shown  to  be 
relatively  small  (a  few  percent),  and  we  have  chosen  to  neglect  these  for  our  simple  models. 
However,  such  effects  can  easily  be  incorporated  in  our  models. 


2.2.  MaRV  Dynamics 

In  this  section,  we  derive  a  simplified  equation  for  the  speed  of  a  MaRV,  and  for  its 
radius  of  curvature  under  a  constant  angle-of-attack  maneuver. 


2.2.1.  MaRV  SPEED 

Above  a  given  height  Hq,  we  assume  that  drag  may  be  neglected,  and  the  MaRV 
travels  at  a  constant  speed  |V^|.  Below  Ho ,  its  speed  decreases  linearly  in  time  until  a 
terminal  (final)  speed  \Vm\  at  or  within  the  keep-out  zone  Hf.0  is  reached.  This  linear 
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deceleration  behavior  is  consistent  with  experimental  tests  and  simulations  [12],  and  is 
shown  in  Figure  2.2  below.  The  vehicle  thus  coasts  at  its  initial  speed  |V^|  until  it  reaches 
the  height  Hq ,  beyond  which  it  decelerates  at  a  rate  of  k^M  feet/sec2  until  the  terminal 
velocity  \V^\  is  reached.  The  time  taken  to  reach  height  Hq  is  simply 


to 


Hq 


XM'ey 


VJK 


(2.4) 


where  is  the  initial  position  of  M  (at  /  =  0)  and  is  its  initial  velocity.  We  assume 
that  Vm  is  constant  in  the  interval  [0,  to]. 


Figure  2.2.  Simplified  speed  profile  of  a  typical  MaRV. 


The  speed  profile  of  M  thus  has  several  segments,  depending  on  whether  the  last  speed 
measurement  is  made  prior  to  reaching  Hq,  or  afterwards. 


'  Wl 

0  <  i  <  to 

|Vj#(t)|=  < 

IVJI 

1  _  (*  -  ro) 

kDM 

,  0  <  To  <  i 

.IVJI 

1  kDM . 

i  To  <  0  <  t 

(2.5) 


2. 2. 1.1.  The  Deceleration  Constant  kjjM 

The  rate  of  deceleration  obviously  depends  upon  the  initial  and  final  positions 

and  velocities  of  the  MaRV,  and  upon  the  height  Ho  where  drag  is  assumed  to  start.  To 
find  kpM)  we  must  first  determine  how  long  M  takes  to  reach  impact  (or  some  other 
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reference  point  inside  the  keep-out  zone),  then  we  must  satisfy  the  initial  and  final  speed 
constraints.  First,  we  solve  for  the  impact  time  tF .  Let  eyo  be  a  unit  vector  in  the  direction 
of  Vjjrf  and  ey  a  unit  vector  along  the  y-axis. 

For  X^L'Cy  >  Hq,  ftDM  =  oo.  Below  Ho, 


X°M-e 


CyO 
VM  y 


—  K\  (i  -  * 


|V"|T„ "  o  “|V"I< 


-TCI 


Applying  the  velocity  constraints  gives:  |V)^|(1  —  tF I^DM )  =  VjjJ,  the  impact  (final)  speed. 


Thus 


Substituting  (2.7)  in  (2.6), 


\K\-\K\ 


evo-ey 


(\Vm\  ~  iyS|)  (|F°|  +  |V£|) 


koM  - 


Hence 

-alV^lgominigo.V^-e,} 

“  (I»SI  -|»J|)  fl»JI  +  IVjSO  ' 

Observe  that  o  is  assumed  to  point  downwards  relative  to  ey,  hence  the  dot  product  in 
the  denominator  is  always  negative  and  kpM  always  positive. 


The  flying  time  from  altitude  Ho  to  final  impact  is  thus 

,F  _ — 2-Hp _ 

0  (lv«l  +  IVjjfl) 


(2.10) 


and  the  total  flying  time  from  to  impact  when  X^»ey  >  Hq  is  TF  =  tq  +  tF ,  where 


Hp-Xj^ey 

eV°-ey 


(2.11) 


When  the  initial  altitude  of  the  MaRV  is  below  Hq,  X^j*ey  <  Hq  and  we  have 

iF  —  ~2  Xjrf^ey 

(\K\  +  I Vm\)  eV°‘ey 


(2.12) 
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In  general,  the  total  flying  time  from  any  starting  position  XQM  is  thus 

TF  =  max{0,  tq}  +  imn{tF ,  tF}  . 


(2.13) 


2.2.2.  MaRV  POSITION 


Starting  from  an  initial  position  X QM,  the  MaRV  position  Xm(<)  at  any  future  time  t 


is 


x m(()  =  X°„  -  Vg,  |i 


t2  —  2  max{0,  ro}< 

2kDM 


max{t  — 


(t 


'  ~  ro,  0}  ] 
-to)  J 


(2.14) 


2.2.3.  MaRV  RADIUS  OF  CURVATURE 


Referring  to  Equations  (2.1)  and  (2.2),  the  instantaneous  radius  of  curvature  of  a  MaRV 


is 


rM(i)  I  = 


2mMexp{*M[*o),e-*— lj 


PoCl(oc,  Mach )mSm 
where  Xm(0  is  the  position  of  the  MaRV  at  time  t. 


(2.15) 


When  banking  or  turning  is  considered  (see  Figures  2.3  and  2.6),  a  banking  or  turning 
angle  <f>\j  must  be  introduced,  and 

2m„exp  { -  l} 

4’m)1  mCL(a,  Mach)„SlV,  ’  (2'16) 

where  the  change  in  height  due  to  turning  through  an  angle  <f>M  in  the  “action  plane,” 
(defined  to  be  the  x-y  plane  of  the  Action  Coordinate  System  (x^,  yfa,  zj^))  is  represented 

by 

A  4>m)  =  ey‘\{Rz&y£Ii&^'>  P'  ^)lrM(<)|) 

x  (sin^M)^  +  (1  -  cos(^M))ejjji)]  *  (2.17) 


where 


,^(7,  (3,  0)  =  R^yR^m-R^iO)  , 


(2.18) 


rotates  the  MaRV  action  coordinates  (xjj,  yjf,  zjfr)  into  the  MaRV  coordinates  (xm,  VM, 
zm)  (the  M-system),  and  <£m|^m=o)  (see  Appendix  A). 
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Figure  2.3. 


The  radius  of  curvature  of  a  vehicle  is  not  only  a  function  of  its  lateral 
acceleration,  but  also  of  its  position  in  the  atmosphere. 


When  the  MaRV  maneuver  is  wholly  within  its  {xm  —  Vm)  plane,  this  rotation  matrix  is 
evidently  the  identity  matrix,  since  then  xj^  =  xm  and  yfo  =  dm- 

It  should  be  clear  from  the  above  discussion  that,  although  we  currently  assume  that 
the  MaRV  maneuver  lies  wholly  within  some  fixed  but  otherwise  arbitrary  plane,  this 
assumption  can  be  eliminated  by  allowing  the  rotation  matrix  R( 7,  /3,  6)  to  vary  with 
time.  Any  MaRV  turn  could  then  be  allowed. 

2.3.  Interceptor  Dynamics 

The  engagement  algorithms  developed  later  in  this  report  allow  the  use  of  two  inter¬ 
ceptors:  an  Overlay  Interceptor  (“overlay”)  and  an  Underlay  Interceptor  (“underlay”). 

Overlays  are  designed  to  intercept  MaRVs  early  during  a  battle  and  at  longer  distances 
from  their  launch  points,  typically  105  feet  to  106  feet,  whereas  underlays  operate  more 
effectively  at  shorter  distances  of  about  104  to  105  feet.  The  former  are  more  sluggish 
and  tend  to  require  more  time  to  reach  their  maximum  velocity.  The  latter  reach  their 
maximum  velocity  considerably  faster,  and  are  designed  for  “close-in”  engagements. 

In  this  section,  we  present  a  summary  of  the  acceleration,  velocity,  and  position  profiles 
of  both  vehicles.  The  complete  derivation  of  these  profiles  is  reported  in  the  appendices. 
We  also  show  the  radius  of  curvature  for  both  vehicles,  as  we  did  for  the  MaRV. 
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2.3.1.  TYPE  I  INTERCEPTORS:  OVERLAYS 

These  are  long-range  two-stage  interceptors  whose  weight  and  thrust  profiles  are  shown 
in  Figures  B.l  and  B.2  and  described  by  Equations  (Bl)  through  (Bll)  in  Appendix  B. 


2. 3. 1.1.  Thrust  Acceleration  of  Type  I  Interceptors 


where 


Using  Equations  (Bl)-(Bll),  from  Appendix  B,  the  acceleration  is 

a/i(0  =  +  Yi(<)a/i2(<)  +  Zi(i)anz(t)  , 

v  u\  _  1  +  sign  (1.61  -t) 

2  ’ 

y^t)  =  si6n(3-83  ~  0^  ^1  +sign(t  -  1.61)^ 


and 


Zi(t)  = 
ain(t)  = 
am(t)  = 


1  +  sign  ( t  —  3.83) 


625,500 


230.5(1  -  0.449t)  +  (m^  -  33.59)  ’ 

_ 114,000 _ 

85.91(1  -  0.159t)  +  (m°n  -  33.59)  ’ 

aJ13(t)  =  0  ,  3.83  <  t  . 


0  <  t  <  1.61  , 

1.61  <  t  <  3.83  , 


(2.19) 

(2.20) 
(2.21) 

(2.22) 

(2.23) 

(2.24) 

(2.25) 


2. 3. 1.2.  Speed  of  Type  I  Interceptors 

Using  Equations  (B16)-(B23),  the  speed  is 

Vn(t)  =  Xi  (t)VIU(i)  +  W)Vm(<)  +  Zi{t)Vm{i)  ,  (2.26) 

where 

VIn(t )  =  — 6043[ln(m®1  +  196.9  -  103.5*)  -  ln(m5,  +  196.9)]  ,  0  <  t  <  1.61  ,  (2.27) 
V112(t)  =  23031n(mjj  +  30.30)  +  60431n(ro?,  +  196.9) 

-  8346  ln(m5!  +  52.32  -  13.66t)  ,  1.61  <  t  <  3.83  ,  (2.28) 


where 


VmW  =  V,°13 


t  -  3.83 
kDIl 


3.83  <  /  <  3.83  +  kDIl 


(2.29) 
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kon  —  10  sec  , 

Vm  =  Vri2(3.83) 

=  23031n(m?i  +  30.30)  +  6043111(771?!  +  196.9)  -  8346 ln(m°i)  .  (2.30) 

A  typical  profile  for  m ?x  =  33.59  is  shown  in  Figure  2.4. 


Figure  2.4.  Speed  profile  of  a  typical  Type  I  interceptor. 


2.3.1. 3.  Position  of  Type  I  Interceptors 

Using  Equations  (B24)-(B30),  from  Appendix  B,  the  distance  traveled  is  (see  Fig¬ 
ure  2.5) 

Xn{t)  =  Xi(t)XIn(t)  +  Y1(t)XIn(t)  +  Zi(t)Xm(t)  ,  (2.31) 

where 

Xlll(i)  =  6043  j  [ln(m?,  +  196.9)  +  l]  t 

+  (^5  +  1-90  -  ^  Mm?!  +  196.9  -  103.52) 

-  ^ - YKs - )  +  196‘9) )  >  0  ^  <  <  1-61  ,  (2.32) 
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XIn(t)  =  X/„(1.61)  +  [2303  hx(m°n  +  30.30)  +  6043  lnfmj,  +  196.9)  +  8346]t 

+  8346  ( +  3-83  -  ln(m7i  +  52.32  -  13-66<) 

\13.66  / 

-  ( 611. Ora® i  +  22,238)  ln^  +  30.30)  +  9729  ln(mj,  +  196.9)  +  13,437  , 

1.61  <  t  <  3.83  ,  (2.33) 


and 

Xm(t)  =  Jfj12(3.83)  +  y;°„ 


3.83  <  t  <  3.83  4-  kjyji 


(2.34) 


Time  (sec) 


Figure 


2.5.  Distance  traveled  for  typical  Type  I  Interceptor. 


2. 3. 1.4.  Radius  of  Curvature 


Similar  to  the  treatment  for  MaRVs  (see  Equations  (2.15)-(2.18),  the  radius  of  curva¬ 
ture  for  both  interceptors  has  the  same  form  and  depends  on  where  the  maneuver  is  made 
(determined  indirectly  by  Xj(Atj)),  and  on  the  turn  angle  (j)j  (see  Figure  2.6). 


MAtj,  <f>i)\  = 


2mJ(AiJ)exp  ^  -  l} 


PoCl(oc ,  Mach )jSj 


where 


Ahj(Atj,  <£/)  =  |r/(A//)|  (1  -  cos  <f)j )  eZM  x  evp  +  (sin  4>i)evo 
and  6 ®  is  the  collision  angle. 


•e 


y  > 


(2.35) 

(2.36) 
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2.3.2.  TYPE  II  INTERCEPTORS:  UNDEF 

These  are  short-range  two-stage  intercepts 
in  Figures  C.l  and  C.2  of  Appendix  C. 


2. 3.2.1.  Thrust  Acceleration  of  Type  II  Interceptors 

Using  Equations  (Cl)-(C15)  of  Appendix  C  the  acceleration  is 

an{t)  =  X2(t)*aI21(t)  +  Y2(t)*aI22(t)  +  Z2(t)*aI23(t)  , 

2  x  106 


where 


aI2l{t)  — 

am(t) 


218.8(1  -  0.714<)  +  m®2  —  15.62  ’ 

105 


109.4(1  -  0.429<)  +  m°I2  -  15.62  ’ 

°J23(0  =  0  ?  2  <t  . 

2.3.2.2.  Speed  of  Type  II  Interceptors 

Using  Equations  (C16)-(C23),  the  speed  is 


0  <  t  <  1 


1  <  t  <  2 


(2.37) 

(2.38) 

(2.39) 

(2.40) 


Vn  (0  =  X2{t)*Vm{t)  +  Y2{t)*VI22(t)  +  Z2(t)*VI23(t)  , 

where 

U/2J (0  =  —12, 803[ln(203.18  +  m°I2  -  156.2f)  -  ln(203.18  +  m°n)],  0  < 

Vi22(t)  =  12,803  ln(203.18  +  m%)  -  1710  ln( 46.98  +  m%) 

-  2131  ln(93.78  +  m°I2  -  46.98t)  ,  0  <  i  <  2, 


Vm(t)  =  Vj°23 


(<~2)\ 
&DI2  J 


2  <  t  <  2  +  kj)i2  , 


(2.41) 
<  1,  (2.42) 


(2.43) 

(2.44) 


where 


k DI2  ~  5  sec  ? 

Vj23  =  12,803  ln(203.18  +  m°I2)  -  10,672  ln(46.98  +  m°I2)  -  2131  ln(m?2)  . 

A  typical  speed  profile  is  shown  in  Figure  2.7. 

2. 3.2. 3.  Position  of  Type  II  Interceptors 

Using  Equations  (C24)-(C31),  the  distance  traveled  equals  (see  Figure  2.8) 

Xj2(t)  =  X2(t)  *  Xj23(t) +  Y2(t)  *  Xj22(t)  +  Z2(t)  *  Xj23(t)  ,  (2.45) 
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Figure  2.7.  Typical  speed  profile  of  a  Type  II  interceptor  with  mass  =  15.62  slugs. 


where 

XI21(t)  =  12, 803 1  (ln(203.18  +  m°I2)  +  1)  +  l)  t  —  ( L3  +  ^|)  M203.18  +  m?2) 

+  (1,3  +  Tw2  ln(203-18  +  m/2  “  156-2  *)}  >  0  <  <  <  1  ,  (2.46) 

Xm (t)  =  XJ21(1)  +  (12,803  ln(203.18  +  m°I2)  -  106.72  ln(46.98  +  m°I2)  +  2131)  t 
-  12,803  ln(203.18  +  m°I2)  -  2131  +  (8541  -  45.44ro?2)  ln(46.98  +  m°I2 ) 

+  (4261  +  45.44m?2  -  2131 1)  ln(93.78  +  m°I2  -  46.9 1)  ,  1  <  t  <  2  ,  (2.47) 


Xm(i)  =  Xm(2)  -  Vf2,  (2  +  ;A- 

V  kDl2 


2  <  t  <  2  +  kf)j2  ,  (2.48) 

where 

I'm  =  3841  ln( 203.18  +  m°l2)  -  1710  ln(46.98  +  m%)  -  2131  ln(m?2). 


2. 3. 2. 4.  Radius  of  Curvature 

The  radius  of  curvature  for  a  Type  II  has  the  same  form  as  that  for  a  Type  I,  with 
possibly  different  values  for  the  payload  mass  mj,  the  lift  coefficient  Ci(a ,  Mach)/,  and 
the  effective  area  Sj. 
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Figure  2.8.  Distance  traveled  for  a  typical  Type  II  interceptor. 


2.4.  Computing  Intercept  Conditions 

The  magnitude  of  miss  distances  depends  very  strongly  on  conditions  which  prevail 
at  the  time  and  location  of  intercept.  In  some  cases,  miss  distances  may  never  have  to 
be  computed  because  of  intrinsic  limitations  in  either  the  MaRV  or  the  interceptor.  For 
instance,  the  MaRV  trajectory  may  not  be  reachable  by  a  given  interceptor,  as  in  the  case 
of  a  MaRV  flying  overhead  at  high  altitude  and  aimed  at  a  distant  target.  Such  cases 
should  be  detected  early  on,  to  avoid  unnecessary  computations  by  the  battle  manager. 

When  a  MaRV  is  reachable,  it  is  important,  to  know — or  at  least  to  estimate — where, 
when,  and  at.  what  angle  of  incidence,  the  interception  will  occur.  If  that  angle  exceeds 
7r/2  radians,  for  instance,  a  tail-chase  condition  arises  where  the  engagement  is  essen¬ 
tially  determined  by  only  the  relative  velocity  of  both  vehicles,  and  relative  maneuvering 
capabilities  become  secondary. 

In  this  section,  we  derive  the  various  intercept  conditions  which  will  be  important  in 
our  miss  distance  estimation. 

2.4.1.  MaRV  TRAJECTORY  REACHABILITY 

If  the  point  on  the  MaRV  trajectory  nearest  to  the  interceptor  launch  site  is  not  reach¬ 
able  by  the  interceptor,  then  surely  no  other  point  on  the  trajectory  will  be.  Consider 
Figure  2.9,  where  XQM  and  are  defined  as  usual,  Xj  is  the  launch  position  of  Intercep¬ 
tor  I,  and  Xq  is  some  intercept  (“collision”)  point  along  the  vector  from  X^j  to  X%j,  the 
MaRV  impact  point.  For  this  report,  we  estimate  by  extending  the  velocity  vector 
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(2.49) 


Vjm  of  M  to  the  (x,  t/)-plaiie.  It  is  easily  shown  that 


V'O 

vF  _  v0  AM’€y  _ 

aM  -  aM - ev°  > 

evo-ey  m 


when  ey  is  a  unit  vector  along  the  y- axis  (when  ev$*ey  =  0,  X^f  obviously  does  not  exist) 


At  the  nearest  point,  X^2 ,  A  Xj  and  evo  are  orthogonal,  i.e .,  their  dot  product  must 


vanish: 


VM 


AXrevo  =  0 

M 


(2.50) 


But  A  Xj  =  Xj/2  -  Xj,  and 


y 


Figure  2.9.  Simplified  interception  geometry  illustrated. 

XC2  =  XM  +  kcevo.  Thus  AXj  =  X°M  +  kcevo  -  xj. 

Hence  (Xj^  +  kcevo  -Xj)»evo  =  0,  or  X0M*evv  +  kc  -  xi*eVo  =  0. 

=  eV&'(X*  ~  x°m)  > 


Thus 


(2.51) 


and 


(2.52) 


Xc/2  =  X°M  +  evo '(Xj  -  X°M)  evo  . 

Observe  that  Xq2  need  not  lie  between  X°M  and  X%j,  as  in  the  case  where  the  interceptor 
is  based  at  location  ( Xj )'.  In  such  cases,  kc  is  negative  and  reachability  is  less  likely. 

To  determine  trajectory  reachability,  we  need  to  solve  for  the  time  t required  by  the 
interceptor  I  to  travel  from  Xj  to  X^2 .  We  thus  need  to  solve 

AX,  =  |*c  2  ~  *°l  =  /  \V,(r)\dr  for  t]g  .  (2.53) 

Jo 

If  there  is  no  solution,  the  MaRV  trajectory — a  fortiori  the  MaRV  itself— is  not  reachable. 
2.4.2.  MaRV  REACHABILITY  AND  INTERCEPT  PREDICTION 

Even  if  the  MaRV  trajectory  is  reachable,  the  MaRV  itself  may  not  be,  because  the 
interceptor  may  be  launched  at  a  time  when  the  MaRV  is  sufficiently  close  to  XfJ2  that 
it  will  have  traveled  beyond  Y^/2  when  the  interceptor  arrives  there.  In  actual  situations, 
of  course,  the  interceptor  guidance  system  will  force  it  to  track  the  MaRV  past  Y^/2,  but 
the  MaRV  may  still  no  longer  be  reachable.  If  it  is,  a  tailchase  ensues  whose  outcome  is 
essentially  determined  by  the  relative  velocity  between  the  two  vehicles. 

Observe  that  the  initial  intercept  estimate  Xc,  as  shown  in  Figure  2.10,  is  not  only 
important  for  reachability  purposes,  but  also  to  establish  initial  conditions  for  problems 
where  interceptor  and  MaRV  maneuvers  are  considered,  as  we  shall  see  in  later  sections. 
Equally  important  is  the  fact  that  the  Xc  is  used  as  an  aimpoint  for  the  interceptor,  which 
is  launched  along  the  vector  A  Xj  from  Y®  to  Xc- 

Consider  thus  a  MaRV  M  traveling  along  a  straight  line  defined  by  V^j  and  by  the  initial 
position  X°M,  and  an  interceptor  I  located  or  based  at  point  Ar°,  as  shown  in  Figure  2.10. 

Let  eyo  and  ej^o  be  unit  vectors  along  the  initial  velocity  V ^  and  the  initial  position 
X°M,  respectively,  d°MC  be  the  distance  from  X°M  to  the  (unknown)  intercept  point  Xc 
and  =  |AY/|  and  A  Xj  =  Xc  —  Xj.  These  quantities  clearly  depend  upon  the  waiting 
time  twa jt  to  launch  the  interceptor  I. 

The  distance  Jmv  =  \X*m  -  X^v\  to  impact  is  simply  Jmv  -  \X%^ey/ev^ey\.  Fur¬ 
thermore, 
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Figure  2.10.  Geometric  conditions  to  derive  the  intercept  point. 


X M  ^Mc( ^wait)eV$  -^c(^wait) 

=  Xl  +  AJV,(W)  ,  (2.54) 

but 

AX/(twait)  =  X%[  +  t^Mc(^wait)eV’°  —  Xj  .  (2.55) 

M 

Hence,  if  we  define  tcw  =  tjc  +  ^wait?  where  tjc  is  the  interceptor  time-to-collision  (total 
flight  time), 

ftlC  ftew 

JQ  \VI(r)\dreAXl=(X0M-X0I)+  I  \VM(r)\dr  evo  ,  (2.56) 

where 


X°M  -  X°I  =  X°MI  =  ( XMIx ,  XMiy,  Xmiz) 


and 


e„.  =  I  ^ 


Vi,.  Vg,y  V ? 


Mz 


M 


n\’  iv&r  n  i 


The  right-hand  side  (RHS)  of  Equation  (2.56)  is  thus 

RHS  =  (x°MU  +  J‘ 


\VM{r)\dr  o  vo 

VMx  >  AMly  + 


l 


V°My 


n\ 


\VM(r)\d.T  , 
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X°MI,+  r*  =(V*(t™)  .  p,(<~)  ,  EM*™))  •  (2-57) 

Jo  KmI  / 

But  the  magnitude  of  RHS  and  the  left-hand  side  (LHS)  must  be  equal.  Hence 
tjc  r  i  1/2 

|V/(r)|dr  =  \Ux(tcw)2  +  Uy(icw )2  +  Uz(tcw)2  J  .  (2.58) 

To  determine  the  time  tjc  elapsed  from  interceptor  launch  to  collision  or  interception,  it 

thus  suffices  to  solve  Equation  (2.58)  for  tic ■ 

'  \ 

j 

To  find  the  root  tic  of  Equation  (2.58),  we  use  algorithm  RTSAFE  [10].  This  algorithm 
requires  the  first  derivatives  of  the  functions  in  Equation  (2.58),  and  these  are  derived  in 
Appendix  D. 


Once  tic — and  hence  tcw — are  known,  d^c  and  <PCI  can  readily  be  found  since 


and 


[tcw 

d°MC  =  /  I  VM(r)\dT 

Jo 

[ tic 

d°ci=  \Vi(r)\dr  . 

Jo 


(2.59) 

(2.60) 


We  can  now  also  compute  the  intercept  angle  0°,  which  is  needed  to  compute  miss  distance 
and  to  determine  the  direction  of  interceptor  launch. 


To  solve  for  0°,  consider  that 


and 


Ac(^wait)  —  ^-M  ^Mci^  wait)eV^ 

A  Xi(twait)  =  Xc(twait)-X0i  . 


(2.61) 

(2.62) 


Hence, 


where  e/^Xj  — 


n  —  AX/(<wait)*eV,°  AXj(t  Wait)*eF° 

C  5  -  lA*fl  "  St,C  |Vj(r)|  dr 

— (-^A I  +  <*Mc(*wait)eVjft  ~  ^’€vm 

~  riwi* 

=  /o‘rc \Vi(T)\dT 

evo,  a  unit  vector  along  the  interceptor  launch  direction. 


(2.63) 
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Chapter  3 

PREDICTING  MISS  DISTANCES 

In  Chapter  1,  we  presented  a  summary  of  the  battle  management  problem  motivating 
this  report.  In  this  chapter,  we  narrow  our  focus  to  the  engagement  of  a  single  MaRV  by 
a  single  interceptor,  as  shown  in  Figure  3.1.  Simple  laws  for  the  behavior  of  both  vehi¬ 
cles  were  derived  in  the  preceding  chapter,  together  with  various  MaRV  reachability  and 
interception  conditions.  First,  we  provide  a  general  overview  of  this  one-on-one  problem. 
Then  we  construct  the  miss  distance  prediction  algorithm  upon  which  the  battle  manager 
is  based,  and  we  discuss  each  of  its  pieces  in  detail. 

In  the  typical  scenario  illustrated  in  Figure  3.1,  a  MaRV  is  headed  towards  some  target 
defended  by  a  collection  of  interceptors.  In  order  to  decouple  our  method  from  major 
tactical  assumptions,  we  do  not  make  the  currently  popular  assumption  that  the  MaRV 
will  execute  its  evasive  or  distracting  maneuvers  at  times  or  altitudes  predetermined  by 
some  internal  computer  program,  although  our  method  accommodates  such  an  assumption 
as  a  special  case.  This  way,  we  also  allow  advanced  MaRV  systems  which  can  estimate  the 
state  of  the  approaching  interceptor  throughout  the  engagement,  either  with  an  onboard 
sensor,  or  from  some  external  source. 

After  acquisition  by  the  Search  and  Acquisition  Sensor  (SAS)  is  complete,  the  Acquisi¬ 
tion  and  Tracking  Radar  (ATR)  tracks  the  MaRV  towards  its  destination.  At  some  optimal 
time  t  =  ti  the  interceptor  is  launched  towards  the  MaRV,  guided  by  the  Guidance  Radar 
(GR)  which  uses  tracking  data  received  from  the  ATR  as  the  engagement  proceeds.  This 
optimal  time  ti  is  the  time  at  which  the  predicted  probability  of  kill  is  a  maximum,  and 
it  is  determined  by  the  predicted  miss  distance,  engagement  altitude,  and  interceptor  war¬ 
head  lethality.  A  capability  to  predict  the  probability  of  kill  achievable  by  an  interceptor 
(the  pk  of  the  interceptor)  is  thus  essential  to  effective  battle  management.  Recall  from 
Chapter  1  that  this  probability  involves  interceptor  lethality  and  interceptor-MaRV  miss 
distance  as  major  parameters: 


Pkq(mocq  =  1)  =  prob{A72ag  >  MDaq}  ,  (3.1) 


where  KRaq 
MDaq 


is  the  kill  radius  of  an  interceptor  of  type  a  against  threat  Oq, 
is  the  miss  distance  achieved  by  an  interceptor  of  type  a  against  Oq. 
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Search  and 
acquisition  sensor 
(SAS) 


Figure  3.1.  Major  elements  of  a  one-on-one  engagement. 


In  this  report,  we  omit  further  discussion  of  lethality  issues,  and  instead  concentrate 
on  the  miss  distance  MDaq.  In  particular,  the  purpose  of  this  chapter  is  to  develop  an 
algorithm  for  predicting  and  minimizing  MDaq ,  and  we  call  this  algorithm  the  MD  algo¬ 
rithm. 

The  miss  distance  achieved  by  an  interceptor  is  determined  by  the  dynamical  state 
of  both  vehicles  at  the  time  of  engagement — their  position,  velocity,  and  acceleration — 
and  by  their  maneuverability.  Because  atmospheric  density  is  an  exponential  function  of 
altitude,  and  the  state  of  both  vehicles  varies  significantly  with  time,  vehicle  dynamics  and 
maneuverability  are  sensitive  functions  of  when  and  where  an  engagement  takes  place.  This 
affords  the  battle  manager  considerable  control  over  the  outcome  of  an  engagement  since, 
by  controlling  the  launch  time  of  an  interceptor,  the  battle  manager  strongly  influences 
the  time  and  place  of  the  engagement,  hence  the  miss  distance  and  the  pk.  Besides  obvious 
resource  parameters  like  interceptor  quantity  and  type,  launch  waiting  time  (tw)  is  in  fact 
the  most  important  parameter  under  the  battle  manager’s  control. 

Interceptor  performance  against  a  MaRV  is  also  sensitive  to  the  quality  of  guidance 
information  and  commands  received  from  the  guidance  radar  and  from  various  sensors.  Not 
only  are  MaRV  state  estimates  corrupted  by  several  types  of  noise,  but  these  estimates 
invariably  reach  the  interceptor  with  some  delay.  This  constrains  the  interceptor  to  base 
its  decisions  on  partially  erroneous  and  obsolete  data,  thereby  increasing  its  miss  distance. 
The  miss  distance  prediction  algorithm  developed  in  this  report  accounts  for  both  sources 
of  error.  The  effects  of  interceptor  system  delays  are  discussed  in  Section  3.4.2,  and  a 
model  for  uncertainty  and  error  propagation  is  presented  in  Section  3.4.3. 

As  shown  in  Figure  3.2,  the  Miss  Distance  (MD)  algorithm  produces  three  principal 
outputs  in  response  to  an  input  Xm  from  the  MaRV  tracker: 

1.  Opportunity  windows. 

2.  Miss  Distance  for  an  undelayed  launch  ( MD° ). 

3.  Launch  direction  vector  (AX/). 

As  the  launch  waiting  or  delay  time  tw  is  varied,  the  miss  distance  function  MD(tw) 
will  typically  exhibit  more  than  one  minimum  for  a  given  MaRV  state  estimate  X^j. 
Regions  where  such  minima  are  attained  are  called  opportunity  windows ,  and  these  are  the 
most  important  outputs  of  MD.  As  we  shall  see  later  in  this  chapter,  such  windows  are 
not  necessarily  connected  subsets  of  the  waiting  time  scale,  and  may  range  in  size  from 
singletons  to  rather  broad  intervals. 

The  particular  miss  distance  MD°  obtained  when  the  interceptor  is  launched  without 


40 


delay  (tw  =  0)  is  important  in  selecting  interceptors  because  it  affects  the  prioritization  of 
targets,  as  we  have  seen  in  Chapter  1.  We  also  have  a  rather  impatient  battle  manager 
that  prefers  launching  too  early  over  launching  too  late,  and  it  uses  MD°  to  select  the 
exact  launch  time. 


Opportunity  windows 


MaRV  state  X„ 
- M  » 

, 

hrbjb  1  =  1,  2,  3} 

Miss  Distance  (MD) 

MD°  =  MD(tw  =  0) 

Technology  parameters^ 

Algorithm 

Launch  direction  vector  AX, 

Figure  3.2.  For  each  MaRV  state  update  Xm ,  algorithm  MD  estimates  opportunity 
windows,  miss  distances,  and  interceptor  launch  direction. 

The  vector  A Xj  pointing  from  the  interceptor  launch  site  Xj  towards  the  predicted 
intercept  point  Xc  is  called  the  launch  vector  and  is  used  as  a  launch  direction  by  the 
battle  manager. 

To  produce  the  three  outputs  shown  in  Figure  3.2,  the  miss  distance  algorithm  analyzes 
an  engagement  in  terms  of  four  regions,  as  illustrated  by  the  typical  geometry  of  Figure  3.3, 
where  Hq  is  the  idealized  boundary  between  the  exoatmosphere  and  the  endoatmosphere, 
and  Hko  is  the  height  of  the  keep-out  zone: 

1.  The  Extremal  Region  (Ro)- 

2.  The  Exoatmospheric  Region 

3.  The  Intermediate  Region  (i^)- 

4.  The  Tailchase  Region  (R3). 

For  each  new  MaRV  position  update  X ^  and  velocity  update  Vjjj  from  the  ATR,  these 
four  regions  are  re-evaluated  to  update  the  opportunity  windows,  assuming  that  the  MaRV 
travels  along  a  straight  line  determined  by  X^j  and  V^.  This  straight  line  continuation 
is  a  good  assumption  because,  even  though  the  MaRV  will  almost  certainly  maneuver  in 
the  future,  the  battle  manager  has  little  knowledge  about  such  maneuvers.  Assumptions 
about  accelerations  and  other  higher-order  terms  are  considerably  less  robust  since  such 
terms  are  subject  to  more  rapid  and  less  predictable  changes.  Furthermore,  as  launch 
time  approaches,  the  time-to-go  decreases,  and  the  possible  error  between  future  MaRV 
velocities  and  their  most  recent  tracker  update  diminishes  rapidly. 
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The  extremal  region  denotes  combinations  of  geometric  and  dynamic  conditions  where 
the  MaRV  is  beyond  reach  of  the  interceptor,  regardless  of  the  value  of  tw.  Consider  for 
instance  an  interceptor  I  based  at  Vj,  and  a  MaRV  M  traveling  from  to  X ^  along  a 

straight  line  trajectory,  as  shown.  Then  M  is  unreachable  by  J,  regardless  of  the  waiting 
time  tun  if  the  point  X^J*  on  the  trajectory  closest  to  Xj  is  not  reachable  by  /. 


Figure  3.3. 


A  simplified  interception  geometry  used  in  the  computation  of  reacha¬ 
bility  and  interception  regions. 


The  exoatmo spheric  region  is  defined  in  terms  of  engagements  which  take  place  at 
altitudes  above  Ho  (typically  100,000  feet).  At  these  heights,  MaRV  maneuverability  is 
severely  limited,  and  the  miss  distance  is  determined  by  the  reach  of  the  interceptor.  If, 
for  some  waiting  time  tw ,  there  exists  a  point  Xf{tw)  in  that  region  which  is  simultane¬ 
ously  reachable  by  both  vehicles — be.,  Xf{tw)  is  a  potential  interception  point — and  if 
the  velocity  of  the  interceptor  arriving  at  the  interception  point  exceeds  zero,  we  set  the 
miss  distance  to  zero  for  this  value  of  tWy  and  tw  is  said  to  belong  to  the  exoatmospheric 
region.  We  thus  assume  that  M  behaves  ballistically  above  Hq  and  that  /,  if  given  per¬ 
fect  information  about  the  position  and  velocity  of  M,  will  intercept  M  without  error  if 
M  is  reachable  for  tw.  The  velocity  constraint  on  I  at  Xf?  is  necessary  to  eliminate  the 
unreasonable  possibility  that  Xf  be  “parked”  at  Xf  until  M  reaches  it. 

The  tailchase  region  is  determined  by  the  relative  velocities  of  the  MaRV  and  the 
interceptor.  If  the  predicted  intercept  or  collision  point  for  a  given  tw  is  Xj^(tw)^  and  if 
the  dot  product  of  their  velocities  at  Xj^(tw)  is  positive,  we  say  that  tw  belongs  to  the 
tailchase  region.  If  the  MaRV  speed  at  the  collision  point  exceeds  that  of  the  interceptor, 
the  MaRV  is  assumed  to  win  the  tailchase,  and  the  miss  distance  is  set  to  zero.  If  not, 
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it  is  set  to  some  very  large  number.  This  condition  for  a  hit  is  a  bit  stronger  than  that 
typically  found  in  the  literature  (see  [13],  for  instance),  but  it  is  simpler,  and  appropriate 
for  our  problem. 

A  waiting  time  tw  belongs  to  the  intermediate  region  if  it  does  not  belong  to  any 
other  region.  It  is  the  only  region  where  evasive  maneuvering  has  a  significant  effect 
on  miss  distance,  and  it  is  thus  also  the  most  difficult  to  analyze  since  the  miss  distance 
algorithm  must  carefully  account  for  the  relative  maneuvering  capabilities  of  both  vehicles. 
As  we  describe  below,  the  miss  distance  estimates  in  this  region  are  based  on  a  simple  but 
revealing  test  consisting  of  a  single  maximal  evasive  turn  by  the  MaRV,  which  is  countered 
by  a  maximal  turn  from  the  interceptor. 


Summarizing  this  introductory  discussion,  the  three  regions  R\ ,  R2,  and  R3  are  subsets 
of  waiting  times  tw,  and  the  MD  algorithm  produces  three  windows  W\  =  [4,11,  412]  C  R\, 
W2  =  [tw 21,  tw22]  c  R2,  and  W3  =  [4,31,  ^32]  C  R3,  as  shown  in  Figure  3.4.  Each  window 
Wi  C  Ri  is  an  interval  of  tw  values  for  which  the  miss  distance  is  a  minimum  in  region 
Ri,  i  =  1,2,  3.  Since  Ro  is  the  extremal  region  where  the  miss  distance  is  arbitrarily  large, 
it  has  no  window,  and  no  information  about  Rq  is  handed  to  the  battle  manager. 


Miss 

Distance 


Figure  3.4.  Algorithm  MD  produces  three  minimal-miss-distance  windows,  one  for 
each  region. 

3.1.  The  Extremal  Region  Ro  (Module  A) 

A  flowchart  for  Region  Rq  (Module  A)  is  shown  in  Figure  3.5. 
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Figure  3.5.  Region  f?o  is  designed  to  eliminate  cases  where  the  MaRV  is  not  reachable 
for  any  waiting  time  tw  (Module  A). 
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Starting  with  zero  waiting  time  (t®  =  0),  the  point  X on  the  MaRV  trajectory  line 
closest  to  the  interceptor  (/)  launch  location  Xj  is  computed  first  (Al).  If  the  MaRV  (M) 
is  moving  away  from  X* ^  (i.e.,  kc  —  ey^o*  ( Xj  —  X <  0),  and  if  I  cannot  reach  M  if 
launched  without  delay  ( tw  =  <® ),  then  I  cannot  reach  M  for  any  tw  (A2,  A3,  A5).  If  M 
is  moving  towards  the  nearest  point  Xc^(kc  >  0),  but  X^'1  is  not  reachable  by  /,  then  I 
again  cannot  reach  M  for  any  tw  (A2  and  A4). 

If  none  of  the  extremal  conditions  hold,  and  the  intercept  point  Xc(t% )  for  <®  is  above 
Ho  (  A6),  then  further  analysis  is  initiated  for  the  exoat mospheric  Region  (i?i)  module  via 
control  line  Ajyg.  If  not,  keep-out  zone  penetration  is  checked  in  Block  A7.  If  )  falls 
in  the  keep-out  zone  (height  Xc(t® )*ey  <  Afc0),  then  the  interceptor  is  ineffective  for  all  tw 
(A8).  If  not,  there  is  no  Region  Rq ,  and  the  lower  limit  <21  of  Region  2  (R2)  is  set  to  t®  in 
Block  A9  (see  our  later  discussion  of  this  region)  and  control  is  passed  to  the  A3  module 
(Module  D)  through  line  A  ad- 


3.2.  The  Exoat  mospheric  or  Nonmaneuver  Region  R\  (Module  B) 

The  purpose  of  this  region  is  to  determine  whether  there  exists  an  opportunity  window 
where  the  miss  distance  is  near  zero  at  exoatmospheric  altitudes  or  in  regions  where  the 
MaRV  cannot  maneuver  (above  Ho).  Referring  to  Figure  3.6,  the  first  operation  in  this 
module  is  to  determine  whether  M  is  reachable  by  I  if  I  is  launched  without  delay  (iw  =  t® , 
Block  Bl).  If  it  is,  MD(f® )  =  0,  t\\  is  set  to  zero  (t®  =  0),  and  flag  1  is  set  since  <u  is  the 
first  value  of  tw  for  which  window  W3  opens  in  Region  3  (B3,  B4).  If  M  is  not  reachable 
(Bl),  the  miss  distance  is  arbitrarily  large  (B2)  and  the  clock  is  incremented  by  the  fixed 
value  of  r  seconds  (B5). 

If  the  current  value  of  tw  is  larger  than  the  maximum  time  tg  allocated  to  Region  1 
(typically  100  seconds),  the  MaRV  has  passed  overhead  and  the  algorithm  stops,  since  M 
never  reached  below  Ho  (B16,B17).  Hence,  there  are  no  Regions  i?2  or  A3 ,  a  fortiori  no 
windows  W2  or  lf'3 .  If  tw  does  not  exceed  if.  and  there  exists  an  intercept  time  tjc{tw) 
above  Ao,  MD(tw)  is  set  to  zero,  as  discussed  in  the  text  (  B16,  B18,  B10,  B9).  If  flag  1 
had  previously  been  set,  the  clock  time  tw  is  again  incremented  by  r  (B8,  B5).  If  not, 
Window  1  opens  at  tw  =  tjj,  and  F\  is  set  (B7,  B6). 

If  there  does  not  exist  an  interception  at  tw  (B10),  MD(tw)  is  set  to  104  (Bll).  If  the 
window  W\  was  previously  closed  (F2  set),  the  waiting  clock  tw  is  incremented  by  r  (B12, 
B5).  If  not.,  but  Wi  was  previously  opened  (F]  set),  we  have  reached  the  end  of  Wy,  ty2  is 
set  to  tw  (B13,  B14),  and  W\  is  now  closed  (B15). 
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Figure  3,6.  Flowchart  for  the  exoatmospheric  region  R\  (Module  B). 


When  the  MaRV  falls  below  H0  (B18)  and  window  W\  was  never  opened  (B19),  there 
is  no  window  in  Region  1  (B20).  If  it  is  still  open,  we  have  reached  the  end  of  Ri,  hence 
the  end  of  W\,  and  the  beginning  t2 1  of  window  W2  (B21).  We  then  analyze  region  R3  by 
passing  control  to  Module  D  via  control  line  Abd- 


3.3.  The  Tailchase  Region  R$  (Module  D) 

Recall  that  we  have  a  tailchase  condition  whenever  V^'(Xc(tw)  -  Xj)  >  0,  i.e.,  any 
time  the  interceptor  velocity  has  a  strictly  positive  component  along  the  MaRV  velocity. 

Module  D  receives  two  inputs,  one  ( Aad )  from  Module  A  and  another  (Abd)  fr°m 
Module  B,  as  shown  in  the  flowchart  of  Figure  3.7.  It  starts  by  resetting  the  flags  previously 
used  in  other  modules  (Dl).  If  the  MaRV  has  reached  the  keep-out  zone  (D2)  prior  to 
being  engaged  in  a  tailchase  (D2),  then  the  current  value  of  iw  equals  the  closing  time  (<22) 
of  window  W2  (D3),  there  is  no  Region  3  (D4),  and  the  waiting  clock  is  set  to  the  starting 
time  <21  for  W2  in  Region  2.  If  there  is  no  W2,  we  are  done  (D6,  D7).  If  not,  we  explore 
R2  with  Module  3  (D26). 

If  M  is  still  outside  the  keep-out  zone  (D2),  but  we  do  not  yet  have  a  tailchase  (D8), 
we  need  to  explore  future  conditions.  To  accomplish  this  efficiently,  an  accelerated  waiting 
clock  is  used  in  which  the  increment  value  is  lOr,  and  not  r  as  before  (D9).  When  the 
first  tailchase  time  is  reached  on  that  “rough”  clock  (D8),  the  waiting  time  tw  is  reset  to 
tw  -  10r,  and  the  start  t22  of  Region  2  is  set  to  tw  (DIO,  Dll).  Using  this  rough  clock  is 
obviously  ten  times  faster  than  using  the  finer  clock  increment  r. 

If  the  speed  of  the  interceptor  exceeds  that  of  the  MaRV  in  the  tailchase  scenario, 
MD(tw)  =  0  (D12,  D13).  If  the  window  W3  was  not  yet  opened  (D14),  it  is  now  opened 
(D16)  and  its  opening  time  <31  equals  tw  (D15).  If  W3  was  previously  opened,  the  clock  is 
incremented  by  r  (D17). 

If  M  is  now  faster  than  I  (D12)  and  if  W3  is  currently  open  (D18),  it  is  now  closed 
(D19)  and  control  is  passed  to  Module  3  to  explore  Region  R2  after  resetting  F\  (D20,  D6, 
D21). 

If,  after  incrementing  the  clock  (D17),  the  MaRV  is  still  above  the  keep-out  altitude 
Hko  (D22),  a  new  relative  velocity  check  is  made  (D12)  and  that  loop  is  re-entered.  If  M 
has  fallen  below  Hko,  and  if  F\  was  never  set,  there  is  no  W3,  and  control  is  transferred  to 
Module  3  (D22,  D23,  D25,  D20,  D21,  D6,  D26).  If  Fi  had  befen  set,  the  amount  of  waiting 
time  tw  equals  the  closing  time  <32  of  Wz  (D24)  and  Module  3  takes  over,  as  before. 
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Figure  3.7.  Flowchart,  for  the  tailchase  region  R3  (Module  D). 
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3.4.  The  Intermediate  Region  R 2  (Module  C) 

In  this  section,  we  discuss  the  most  complex  engagement  region,  region  Ri.  Because 
the  algorithms  dedicated  to  this  region  are  based  in  part  on  heuristic  arguments  which  are 
not  obvious,  we  have  devoted  a  special  introductory  section  (Section  3.4.1)  to  presenting 
and  defending  these  arguments.  Then,  in  Section  3.4.2,  we  develop  the  mathematical 
machinery  required  to  define  miss  distance  in  Region  2.  Finally,  in  Section  3.4.3  we  present 
and  discuss  the  flowchart  for  Module  C,  as  we  did  in  previous  sections  for  the  other  regions. 

3.4.1.  INTRODUCTION 

Two  principal  characteristics  distinguish  the  intermediate  region  (R2)  from  the  other 
regions.  First,  engagements  in  this  region  take  place  endoatmospherically  (below  Ho), 
where  the  atmosphere  is  sufficiently  dense  that  drag  and  lift  become  significant,  and  evasive 
maneuvering  is  the  principal  method  for  avoiding  interception.  Secondly,  and  by  definition, 
tailchasing  does  not  arise  in  this  region,  and  the  closure  rate  between  MaRV  and  interceptor 
is  thus  positive. 

Since  maneuvering  capabilities  are  such  strong  determinants  in  this  region,  and  little 
is  known  about  the  future  maneuvering  tactics  of  a  MaRV,  how  then  can  miss  distances  be 
accurately  predicted  in  this  region?  This  question  has  obviously  long  been  of  significant 
interest  to  many  analysts  and  interceptor  designers,  and  several  analyses  and  methods 
[13-19]  have  been  reported. 

In  the  context  of  battle  management,  one  popular  method,  summarized  in  [19],  is  based 
on  a  combinatorial  approach  where,  on  the  basis  of  intelligence  data  and  various  technolog¬ 
ical  constraints,  representative  MaRV  maneuvers  and  trajectories  are  derived.  Interceptor 
capabilities  against  such  trajectories  are  then  evaluated  using  computer  simulation  meth¬ 
ods.  In  spite  of  its  popularity,  disadvantages  of  this  method  appear  obvious,  particularly 
for  battle  management  purposes.  First,  there  is  an  uncountable  infinity  of  trajectories, 
and  to  evaluate,  design,  or  deploy  an  interceptor  on  the  basis  of  just  a  few  of  these  is 
very  risky.  Even  if  these  few  contain  some  “optimal”  trajectory,  this  approach  encourages 
gaming  by  the  adversary,  and  may  result  in  his  selecting  a  trajectory  which  is  suboptimal 
by  the  assumed  criterion,  but  optimal  by  another — equally  rational — criterion.  Second, 
this  method  does  not  sufficiently  account  for  real-time  information  that  is  available  from 
sensors  as  the  engagement  evolves.  Reliable  state  information  about  a  MaRV’s  progress 
towards  its  target  for  instance  will  significantly  reduce  its  maneuvering  choices,  particu¬ 
larly  as  it  nears  the  target.  Third,  recall  that  the  principal  reason  why  the  battle  manager 
needs  miss  distance  information  is  to  effectively  allocate  interceptor  type  and  quantity, 
launch  time,  and  launch  direction.  Traditional  miss  distance  prediction  methods  are  of 
little  assistance  in  making  such  decisions. 
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Another  approach  (see  Gavel  [20],  for  instance)  employs  optimal  control  methods  to 
obtain  some  “optimal”  trajectory  which  can  be  used  in  simulations  to  derive  miss  distance 
estimates.  While  this  method  is  still  under  investigation,  it  appears  to  have  limitations 
relating  to  the  choice  of  a  criterion  and  to  the  uniqueness  of  a  solution,  given  that  engage¬ 
ment  dynamics  are  highly  nonlinear  in  this  region.  Furthermore,  expressing  the  problem 
as  a  multidimensional  optimization  problem  has  severe  real-time  computational  limita¬ 
tions,  particularly  for  engagement  scenarios  where  the  MaRV  does  not  in  fact  follow  a 
precomputed  optimal  trajectory,  and  real-time  adaptation  must  be  considered. 

Another  method  is  popular  in  the  guidance  community,  one  which  uses  linear  dynamical 
models  of  both  vehicles  to  estimate  miss  distances,  as  exemplified  by  Reference  [14]. 

In  contrast  to  these  methods,  our  approach  is  based  on  an  evasion-and-capture  test 
whose  outcome  is  determined  by  the  maximal  technological  capabilities  attributed  to  both 
vehicles.  Whereas  other  methods  rely  on  detailed  tactical  assumptions  and  on  specific 
guidance  algorithms  to  estimate  miss  distances,  and  could  thus  be  classified  as  “scenario- 
based,”  we  have  intentionally  avoided  such  assumptions  in  developing  our  test,  and  our 
approach  should  probably  be  classified  as  “technology-based.”  Its  basic  premise  is  rather 
simple,  and  can  be  characterized  by  a  single  question: 


Assuming  a  “reasonably  rational  and  informed  MaRV-Interceptor  pair,”  what,  is 
the  “best”  that  the  interceptor  can  do  to  capture  the  MaRV,  given  the  current 
state  and  environment  of  both? 


To  see  what  we  mean  by  “reasonably  rational  and  informed”  and  by  “best,”  refer  to 
Figure  3.8,  where  we  show  a  typical  scenario  in  which  a  MaRV  attempts  to  evade  an  in¬ 
terceptor  approaching  at  an  angle  0C  with  a  single  right-hand  turn,  which  the  interceptor 
counters  with  a  corresponding  left-hand  turn.  We  say  that  the  MaRV-Interceptor  pair  is 
reasonably  informed  if  they  both  possess  some  state  information  about  each  other  through¬ 
out  the  engagement.  In  the  case  of  the  MaRV,  this  assumption  deserves  some  discussion 
because  the  MaRV  is  conventionally  assumed  to  operate  in  an  open-loop  mode,  unaware 
of  the  interceptor.  As  we  shall  see,  this  is  not  important  in  our  minimax  launch  logic  be¬ 
cause  an  interceptor  is  launched  at  the  precise  time  at  which  the  maximum  miss  distance 
is  a  minimum.  The  location  of  this  minimum  is  determined  by  atmospheric  and  dynamic 
conditions  unrelated  to  the  information  state  of  the  MaRV. 

Although  the  location  of  the  minimum  miss  distance  on  the  launch  time  scale  is  in¬ 
sensitive  to  the  MaRV  information  state,  its  magnitude  generally  is  not.  This  may  lead 
to  considerably  overestimating  the  size  of  the  minimum  miss  distance,  particularly  if  the 
MaRV  is  known  not  to  maneuver  at  all.  Such  an  error  would  lead  the  interceptor  man- 


50 


ager  to  underestimate  the  probability  of  kill,  and  might  cause  an  excessive  expenditure 
of  interceptors  in  special  cases  so  favorable  to  the  defense.  But  such  MaRV  disabilities 
are  virtually  impossible  to  predict  in  practice,  hence  have  little  operational  value  from  a 
defense  allocation  point  of  view.  Whether  the  MaRV  has  interceptor  state  information  or 
whether  it  follows  a  preprogrammed  trajectory,  the  defense  cannot  predict  its  maneuver, 
and  must  assume  that  it  will  maneuver  to  maximize  its  distance  from  the  interceptor.  It 
is  this  maximum  distance  which  is  minimized  in  the  MD  algorithm. 


Figure  3.8.  Illustrating  a  single-turn  evasive  maneuver. 
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Delays  in  interceptor  response  can  be  major  contributors  to  miss  distance.  Not  only 
is  MaRV  state  information  delayed  when  it  arrives  at  the  interceptor,  but  the  guidance 
and  control  subsystem  of  the  interceptor  cannot  respond  instantaneously.  In  Section  3.4.2 
below,  we  discuss  how  the  MD  algorithm  accounts  for  the  causes  and  effects  of  such  delays. 
State  information  received  by  the  interceptor  can  also  be  severely  distorted  by  various  error 
mechanisms  in  the  propagation  and  processing  of  information.  To  account  for  the  effects  of 
random  errors  and  information  distortion,  we  develop  a  sensitivity  and  error  propagation 
model  in  Section  3.4. 

We  say  that  a  MaRV  and  an  interceptor  are  reasonably  rationale  neither  maneuvers  too 
early  or  too  late,  based  on  their  knowledge  of  each  other’s  state.  Referring  to  Figure  3.8, 
a  MaRV  following  trajectory  T\  has  turned  too  early,  since  it  gives  the  interceptor  an 
opportunity  to  “goal-tend,”  a  maneuver  providing  significant  leverage  to  the  interceptor. 
If  M  delays  its  evasive  maneuver  too  much,  and  takes  trajectory  T3,  another  advantage 
is  given  the  interceptor  since  the  latter  need  not  turn  as  much  or  as  quickly  to  intercept 
M .  A  dramatic — and  admittedly  extreme — example  of  this  advantage  is  where  M  delays 
its  turn  so  long  that  it  would  forfeit  any  opportunity  to  maneuver  before  intercept,  which 
would  take  place  at  Xc.  Normally,  this  would  arise  only  if  M  behaved  ballistically,  as  if 
it  were  an  ordinary  RV.  From  a  trajectory  point  of  view,  a  similar  advantage  is  given  the 
MaRV  if  the  interceptor  turns  too  late.  The  MaRV  then  escapes  along  trajectory  T3. 

So  the  only  scenario  which  is  equitable  to  both  the  MaRV  and  the  interceptor  is  the 
second  trajectory  (T2),  which  touches  the  interceptor  trajectory  at  one  point  Q.  The 
location  of  both  vehicles  along  their  respective  locus  at  any  given  point  in  time  obviously 
depends  on  their  velocity  and  the  system  time  delays  experienced  by  the  interceptor.  Using 
the  above  assumptions,  the  miss  distance  is  thus  determined  only  by  the  system  time  delay 
and  by  the  relative  velocity  and  acceleration  capabilities  of  the  vehicles,  and  not  by  tactics 
or  strategy,  and  that  is  why  we  describe  our  approach  as  technology-based. 

An  equally  compelling  reason  why  both  the  MaRV  locus  and  the  interceptor  locus  must 
meet  at  only  one  point  follows  from  the  requirements  for  mathematical  consistency.  Clearly, 
if  both  vehicles  have  exactly  the  same  capability  and  errors,  and  delays  are  neglected,  the 
miss  distance  should  be  zero  since  we  do  not  allow  tactical  advantages.  Futhermore,  if 
the  velocity  or  acceleration  capabilities  of  the  interceptor  are  increased,  the  miss  distance 
could  not  possibly  get  worse,  and  should  remain  at  zero.  But,  if  the  two  trajectories  (loci) 
do  not  remain  in  contact  as  interceptor  capabilities  are  increased,  this  is  impossible.  Since 
the  late-maneuver  (T3)  was  disallowed  earlier,  the  two  loci  must  meet  at  one  and  only  one 
point. 

Having  developed  and  defended  our  fairness  doctrine  for  a  single-turn,  we  now  proceed 
to  show  that  the  outcome  of  an  entire  engagement — and  hence  the  miss  distance  between 
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the  MaRV  and  the  interceptor — may  be  predicted  from  this  simplified  engagement. 

Regardless  of  its  complexity,  a  trajectory  is  a  sequence  of  turns,  potentially  infinite. 
It  is  well  known,  however,  that  considerable  energy  and  momentum  is  lost  during  a  ma¬ 
neuver,  and  that  long  and  complex  maneuvers  are  not  effective.  What  current  technology 
encourages  is  a  trajectory  consisting  of  one  or  two  carefully  timed  turns,  the  second  of 
which  typically  takes  place  at  reduced  velocities  and  altitudes  where  the  interceptor  usu¬ 
ally  has  a  speed  or  maneuvering  advantage.  Considering  that  the  battle  manager  employs  a 
two-layer  shoot-look-shoot  strategy  where,  if  a  MaRV  is  missed  during  a  first  engagement 
with  an  overlay  interceptor,  an  underlay  interceptor  is  launched  to  engage  it  during  its 
second  maneuver,  this  maneuvering  advantage  further  reduces  our  concern  for  the  MaRV’s 
second  maneuver.  The  relative  capabilities  of  both  vehicles  can  thus  be  estimated  using 
very  simple  maneuvers,  sometimes  a  single  turn  [16]. 

Recall,  however,  that  any  engagement  outcome  is  strongly  influenced  by  where — and 
thus  when — the  engagement  takes  place.  A  single  test  at  a  fixed  altitude  would  not  prop¬ 
erly  represent  all  reasonable  possibilities,  and  would  be  comparable  to  the  combinatorial 
methods  which  we  criticized  earlier.  Instead,  our  algorithm  sweeps  all  altitudes  via  the 
waiting  variable  tw ,  to  determine  when  the  predicted  miss  distance  is  a  minimum.  At  this 
minimum,  called  t*w,  the  interceptor  is  said  to  do  best.  We  thus  literally  examine  an  infinite 
set  of  trajectories  and  engagements  before  settling  on  a  miss  distance.  We  claim  that  this 
infinite  set  adequately  represents  all  reasonable  MaRV  maneuvers. 

3.4.2.  MATHEMATICAL  ANALYSIS  OF  REGION  2 

Now  we  develop  the  additional  machinery  required  to  derive  miss  distance  predictions 
for  Region  i?2-  Referring  to  Figure  3.9,  our  objective  is  to  find  an  expression  for  the  min¬ 
imum  distance  between  the  MaRV  and  the  interceptor,  as  both  are  negotiating  maximal 
turns.  Roughly  speaking,  our  algorithm  involves  three  major  steps.  In  the  first  step,  it 
computes  the  interception  point  Xc  based  on  the  continuation  of  the  MaRV  velocity 
The  derivation  of  Xc  was  discussed  in  detail  in  Section  2.4.2.  The  second  step  is  to  con¬ 
struct  the  MaRV  and  interceptor  loci.  Although  we  show  the  two  loci  touching  each  other 
in  Figure  3.9,  this  iterative  construction  process  starts  with  a  separation  A,  as  illustrated 
in  Figures  2.6  and  3.10,  and  is  designed  to  reduce  A  to  near  zero  in  a  minimal  number 
of  steps.  Because  of  the  severe  atmospheric  and  dynamical  nonlinearities  in  our  vehicle 
model,  transcendental  equations  must  be  solved.  Since  such  equations  have  no  satisfactory 
closed-form  approximations,  we  are  forced  to  employ  numerical  optimization  methods  to 
solve  for  the  unique  point  where  both  trajectories  meet,  when  their  separation  A  has  been 
reduced  to  zero.  Given  that  the  starting  point  in  any  nonlinear  optimization  problem 
is  crucial  to  the  accuracy  and  efficient  convergence  of  the  optimization,  a  considerable 
amount  of  effort,  was  devoted  to  choosing  this  initial  condition.  By  neglecting  atmospheric 
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nonlinearities  for  this  initial  condition,  we  were  able  to  exploit  the  fact  that  both  vehicles 
travel  along  the  perimeter  of  a  circle,  and  obtained  a  closed-form  solution  which  is  exact 
when  atmospheric  effects  are  neglected. 

More  formally,  consider  Figure  3.11,  where  all  the  initialization  parameters  are  illus¬ 
trated.  For  this  approximation,  the  loci  of  both  vehicles  are  circles  which  are  assumed  to 
touch  only  at  one  point.  The  distance  Afj,  0® )  between  M  and  I  at  the  time 


Figure  3.9.  Miss  distance  geometry. 
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A i°M  =  At0  +  tw  when  both  are  initiating  their  maneuvers  (6  =  0  for  this  initialization) 
and  when  the  collision  angle  is  6°,  is  then 

dMi(tw ,  At0!,  0°)  =  2  [|rM(At^)||r/(At?)|]1/2  cos  6°  +r/(Af^)  sin(0°)  .  (3.2) 


Figure  3.10.  Geometric  relationship  between  optimization  variables  showing  the  ini¬ 
tial  offset  A. 


Obtained  by  purely  geometric  considerations,  this  distance  must  equal  the  distance 
R°mi  —  A Rmi  separating  the  two  vehicles  after  M  has  traveled  for  A t°M  =  tw  +  A t® 
seconds  and  I  has  traveled  for  At®  seconds,  where  R^j  is  defined  in  Figure  3.11,  and 

/•A  t^+tw  »Atj° 

ARMi(tw,  At°i)  =  /  \Vm(t)\cIt  +  /  \Vj(r)\dT  .  (3.3) 

Jo  Jo 
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Figure  3.11.  Parameters  and  geometry  for  the  initialization  process. 
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But  At®  is  not  known,  hence  must  be  solved  for..  So  the  initialization  of  the  control  variable 
A tj  may  be  reduced  to  solving  the  equation 

Fo(tw,  At?,  e°c)  =  dMi(tw ,  At/,  6°c)  -  R°MI(tw )  +  a RMI(tw,  At?)  =  0  (3.4) 

for  At/.  The  solution  to  this  equation  is  found  by  using  the  root  solving  procedure  RTSAFE 
([10]  and  Appendix  F)  and  is  defined  as  At?.  • 

Starting  with  this  initial  estimate,  all  that  remains  is  to  adjust  both  loci  to  account 
for  atmospheric  effects,  and  this  is  done  with  Algorithm  DBRENT  ([10]  and  Appendix  E), 
and  with  a  modified  form  of  the  Newton-Raphson  approach.  Realizing  that  the  principal 
independent  variable  to  slide  the  two  loci  towards  or  away  from  each  other  is  the  time-to-go 
A  tj  (assuming  the  waiting  time  tw  is  fixed),  the  algorithm  computes  the  partial  derivative 
dA(tw)/dAti,  where  A(tw)  is  the  distance  between  both  loci.  Then,  after  computing  A(tw) 
using  DBRENT,  a  correction  value  d(Ati(tw))  =  A(tw)[dA(tw)/dAti]  1  for  Ati(iw)  is 
used  to  adjust  the  time-to-go  accordingly,  thereby  sliding  the  loci  to  reduce  A. 

This  iterative  process  continues  until  A  falls  below  an  acceptance  threshold  e,  which 
was  arbitrarily  set  to  10  feet.  Because  all  the  required  derivative  information  is  available 
in  closed  analytical  form  (see  Appendix  E),  this  process  converges  very  rapidly;  in  a  single 
step  for  the  vast  majority  of  cases  tested.  In  some  rare  instances,  several  steps  are  needed, 
and  we  have  set  a  limit  of  5  on  the  number  of  steps  allowed.  The  penalty  resulting  from 
this  limit  has  never  exceeded  10.0  feet  of  error  in  our  miss  distance  tests. 

The  two  loci  now  joined,  the  third  and  final  step  in  the  procedure  is  to  compute  the 
miss  distance.  Referring  again  to  Figure  3.9,  the  MaRV  initiates  its  turn  at  the  “banking 
time”  f/j,  and  the  interceptor  at  a  delayed  time  tjg  +  8,  where  8  represents  the  various 
interceptor  time  delays.  In  this  example  the  MaRV  reaches  Q*  ahead  of  the  interceptor, 
at  time  <2  =  The  interceptor  reaches  Q*  at  time  t$  =  tj,  when  the  MaRV  is  ahead  at 
position  Xj^(tj).  The  miss  distance  estimate  is  now  defined  as 

f  {\X M(l  +  8)  -  A"/(f)|}  ,  t*M  <  t*j  +  8 

MD(tw)  =  {  ,  (3.5) 

U,  +  S 

where: 

tg  is  the  time  at  which  the  MaRV  initiates  its  turn, 
t*M  is  the  time  taken  by  M  to  reach  Q* , 
tj  is  the  time  taken  by  I  to  reach  Q *, 

8  is  the  interceptor  system  time  delay. 

To  provide  a  helpful  comment,  the  miss  distance  should  be  zero  if  the  interceptor  reaches 
Q*  before  the  MaRV  does  since  the  interceptor  can  then  turn  towards  the  MaRV  and 
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intercept  it  before  it  reaches  Q*.  If  the  MaRV  gets  to  Q*  first.,  the  validity  of  the  first 
expression  is  rather  obvious,  since  MD(tw)  can  only  get  larger  beyond  Q*.  As  a  good 
approximation  to  MD{iw)  in  this  condition,  and  to  save  computation  time,  Algorithm  MD 
actually  computes  the  simpler  expression 

MD(tw)  =  \Q*  -XM(t})\  ,  (3.6) 

where  X^[(t^)  is  the  position  of  M  at  the  precise  time  t*j  when  I  reaches  Q* . 

For  a  more  rigorous  derivation  of  the  distance  A (tw)  between  both  trajectories,  refer 
again  to  Figure  3.10,  and  consider  some  arbitrary  point  Q  on  M’s  trajectory  or  locus. 
Then 

A(tw)  =  min  {[\r'j (<f>i I  -  \rl{<t>l{<t>M))  I]}  ,  (3.7) 

4>M 

where  <f>i((/)M)  is  the  angle  swept  by  I  as  it  travels  from  the  start  of  its  turn  at  Xj(Atj) 
to  the  point  where  its  extended  radius  of  curvature  passes  through  Q ,  as  shown.  Clearly 
fili&M )  is  indeed  determined  by  <f>M  only.  We  define  as  <f>*M  the  MaRV  angle  <f>M  at  which 
this  minimum  occurs.  Observe  also  that  this  minimum  could  be  negative,  in  which  case 
the  loci  must  be  pulled  away  from  each  other. 

Using  now  the  MaRV  coordinate  system  M  =  VMi  zm)i  and  standard  vector 

analysis  methods, 

Q  =  XM(AtM)  -  \rM(AtM)\eyM  +  tm{4>m)  (3.8) 

and 

Q  =  XM(AtM)  +  dMc{tw ,  A tj)eXM  -  da(tw,  A ti)eAXj  +  ri{Atj)eZM  X  eAXr 

+  r'l{<t>M)  ,  (3.9) 

where: 

A tj  is  the  interceptor  flight  time, 

tw  is  the  launch  waiting  time, 

At m  =  tw  +  Atj  is  the  MaRV  flight  time, 

AV/  is  the  interceptor  launch  direction  vector, 

is  the  angle  swept  by  M  as  it  travels  from  Xm( At^) 
to  Q. 


tm(4>m)  =  Vm{4>m)\ (cos  (f>M)eyM  +  \tm(<!>m)\ (sin  4>m)cxm  ,  (3.10) 

eAXj  =  —(cos 0c)exM  -  (sin BQc)eyM  =  evo  ,  (3.11) 


and 
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\rj(Ati)\eZM  x  evo  =  |rj(Aij)|  (sin  0%)eXM  -  |r/(  A</)|  (cos 0°c)eVM  ,  (3.12) 

where: 

0®  is  the  incidence  (collision)  angle, 

tm{4>m)  is  the  curvature  radius  of  M  as  a  function  of  its  angular 
displacement  4>Mt 

rj(Atj)  is  the  curvature  radius  of  I  as  a  function  of  its  flight  time  At], 
TM(AtM)  is  the  curvature  radius  of  M  as  a  function  of  its  flight  time  Atjif. 


Since 


and 


eZM  X  eyp  =  (sin0°)eXM  -  (cos0°)ej,M  ,  (3.13) 

/•A 

dMc{tw,  A tj)  =  \Xm(A tjvf)  ~  Xc(tw)\  =  \X°M  —  Xc(tu,)|  —  /  | Vm(t)\(1t  (3.14) 

Jo 

/•Atfcf 

dci(tw,  &ti)  =  \Xc(tw)  -  X/(Af/)|  =  |Jfj  -  Xc(tw)\  -  /  \V,(T)\dr  .  (3.15) 

Jo 


Hence 


Q  =  Xm{AIm)  -  \rM(AtM)\eyM  +  \tm(<I>m)\ (cos (^mVvm  +  \tm{4>m)\  (sin  <f>M)eXM 
=  Xm{AIm)  +  \r  (sin  4>M)eXM 

+  [|rM(0Af)|  COs(4>M)  ~  \rM(AtM)\]eyM  (3.16) 

Therefore: 

r\{<i>M)  =  -VM{^M)\^yM  +  \rM{<l>M)\(™s<f>M)eyM  +  \rM{<t>M)\ (sin  </>M)eXM 

wait?  ^ i)€xm  ~~  d'Cj(t wait?  A/j)  (cos  9C  )eEM  +  \r At j)  \  (cos  Oc)eyM 
~dci{t wait?  Aij)  (sin  0C )e3/Af  —  |r/(A^/)|  (sin^c 

=-4^xAf  +  BeyM  ,  (3-17) 

where 

A  =  (shi  4>m )  ^Mc(^wait  7  Atj)  ~  ^C/(^wait?  (cos^c  ) 

-  |rj(A</)|(an<")  .  (3.18) 

and 
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B  =  (cos  (f>M )  -  \rM{AtM)\  -  dc/(<wait,  At j)  (sin  0 {) 

+  MA</)|(costfJ)  . 

To  obtain  rj(4>M)i  it  now  suffices  to  compute  4>j{4>m)  and  rI 


To  find  <f>j, 


cos  (<f>j) 


~rJ*(e2A f  x  fAjfj)  _  rJ‘e*M  x  eVJ° 


(3.19) 


_  rj(eZM>  CS/An  C^m)*  (e^M  X  (COS^c)  eXM  “I"  (sin  0c)eVM ) 


_  rl(egAM  eyjtf>  e^Af)*  (COS^c)  eVAf  (Sin  ^c)  el 


_  (^eZAf  +  (cos 0|? )  eyjtf  -  (sin0°)  ea 


-A  (sin  0°)  +  5  (cos 


Hence, 


cos(^j)  = 


—  A  (sin  0°)  +  B  (cos  0°) 


To  find  rj  (4>j  (4>m)),  we  now  simply  use  Equations  (2.27)  and  (2.28). 


(3.20) 


(3.21) 


To  compute  the  miss  distance,  we  need  the  time  tj  at  which  the  interceptor  reaches 
Q* *  But  the  angle  <f>j  swept  by  I  to  reach  Q*  was  found  by  numerical  methods  earlier 
(Equations  (3.2)-(3.4),  and  with  Equation  (3.21).  We  can  thus  find  t*j  by  solving 

•  (3-22) 


To  avoid  a  difficult  nonlinear  integration,  we  approximate  r/(r)  by  [rj(tg)  x  ?*/(</>})] ■ 
This  is  a  very  good  approximation  since  rj  does  not  vary  much  during  a  sharp  turn, 
particularly  at  low  altitudes  where  rj  is  small  and  thus  the  loss  of  altitude  during  a  bank 
will  also  be  small.  This  function  V/(r)  and  its  integral  are  derived  in  Appendix  B. 

Continuing,  if  the  interceptor  reaches  Q*  before  the  MaRV  does,  the  miss  distance  is 
zero,  since  it  can  simply  turn  towards  the  MaRV  as  it  approaches  Q* .  Mathematically,  if 


Atu+tRf+i 


VM(r)di 


'A^=At/°  +*-  [rM{AtQM)rM{<f>*M)\ 


<  <j>) ,  then  MD  =  0 


(3.23) 
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If  not,  we  need  to  compute  the  time  11bm  required  by  M  to  reach  Q* ,  by  solving 


4>*m  -  [ 

JA 


AtM+tBM  Vm(t)(1t 


AtM  rM(At%)rM{(f>*M)] 


for  t*BM  and  evaluating 


AiM+tBM 


Vm{t)  dr  . 


A  closed-form  solution  for  t*nM  can  be  obtained  as  follows: 


fAtM+  tBM 
=  /  ^l[l  —  ^2T]  dr  —  k\T 


|  At0M+fBM  klk2T2\ 


(3.24) 


(3.25) 


=  kit*BM - +  2A 


In  standard  quadratic  form: 


t*BM  +  (*4*1*2  _  *iKbM  +  <j>*M  —  0  . 


(3.26) 


Solution: 


(1  -  A  t%k2)  -  J(1  -  A  t%k2f  -  2f^M 


(3.27) 


where 


H0-X°M-ey 

lf  To  “  Fo.e  <  0  ’ 

vm  ev 


and  if  tq  >  0, 


ki  = - - — — — - 7-  and  k2  =  — - 

[rM(Xt0M)rM(<t>M)} 1/2  RDM 


\K\  [X  + 

[r  M  ( A4  )r  m  (  4>*m  )] 


]/2  "2  [kDM  +  To] 


3.4.3.  UNCERTAINTY  ANALYSIS  AND  ERROR  PROPAGATION 


Although  parameter  and  state  vector  errors  can  affect  miss  distance  estimates  in  all 
three  regions,  Regions  1  and  3  are  rather  insensitive  to  such  errors  because  of  their  binary 
nature.  Miss  distances  tend  to  be  either  very  high  or  very  low  in  these  regions,  and  only  the 
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width  of  their  windows  can  be  affected.  Battle  management  decisions  are  rarely  sensitive 
to  this  width. 

In  Region  2,  however,  state  vector  or  parameter  errors  can  be  significant,  and  in  this 
section  we  develop  a  model  with  which  their  effects  can  be  quantitatively  assessed. 

Even  though  miss  distance  prediction  in  Region  2  requires  the  application  of  numerical 
methods  because  the  overall  problem  has  no  closed  form  solution,  once  the  MaRV  and 
interceptor  loci  are  glued  together,  the  problem  becomes  formally  tractable  once  again. 
The  various  derivatives  obtained  in  the  appendices  can  now  be  employed  as  sensitivity 
analysis  and  error  propagation  tools,  if  some  mild  assumptions  are  made  about  the  nature 
of  the  input  errors,  as  symbolized  in  Figure  3.12. 


[(xM>  e)  cov(x  Mi  01 
Input  errors 


Miss  Distance 
(MD) 
Algorithm 


[(MP,  AXt,  tL)  cov(MP,  AX1(  tL)]^ 
Output  errors 


Figure  3.12.  Algorithm  MD  computes  the  sensitivity  of  output  variables  to  errors  and 
uncertainties  in  input  variables. 


Consider  the  MaRV  state  Xm  and  the  parameter  vector  6  corrupted  by  errors  and 
uncertainty  summarized  by  the  covariance  matrix  cov(Xji/ ,  6).  If  we  assume  that  errors 
are  small,  Gaussian,  and  additive,  and  that  MD  is  linear  for  small  input  disturbances,  the 
covariance  of  the  Gaussian  output  random  vector  (MD,  AX/,  f/)  is  simply 

co v(MD,  AX/,  <2)  =  J  cov(Xaf,  0)JT  ,  (3.28) 

where 

_  djMD ,  AX/,  tj) 

d{XM,  0) 

the  Jacobian  of  the  output  with  respect  to  the  input  of  MD.  Since  the  partial  derivatives 
in  J  are  available,  the  computation  of  the  output  error  distribution  is  routine. 

3.4.4.  MISS  DISTANCE  MODULE  FOR  REGION  R2  (MODULE  C) 

Recall  that  a  waiting  time  tw  falls  in  region  2  if  it  falls  in  no  other  region.  In  contrast  to 
the  other  regions  where  miss  distances  are  either  very  small  or  very  large,  miss  distances  in 
R.2  usually  assume  a  continuum  of  values,  ranging  typically  from  a  few  feet  to  several  thou¬ 
sand  feet.  Because  this  region  is  computationally  complex,  we  do  not  actually  compute 
the  entire  window  W2,  but  only  find  the  optimal  waiting  time  t*w  and  the  minimum  miss 
distances  MD*  at  that  optimal  time.  This  approach  is  perfect  if  the  minimum  is  unique, 
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since  the  window  has  zero  width  and  t fully  specifies  this  singleton- window.  Sometimes, 
a  flat  spot  in  the  MD  function  does  exist,  but  this  is  not  important  in  practice  because 
such  spots  are  typically  very  narrow.  Exceptions  usually  arise  in  those  cases  where  the 
interceptor  totally  dominates  the  MaRV,  but  in  such  cases  the  battle  manager  will  have 
ample  opportunities  to  neutralize  the  MaRV  near  that  flat  spot,  or  in  the  other  opportu¬ 
nity  windows  which  will  likely  arise.  A  simplified  flowchart  of  Module  C  is  illustrated  in 
Figure  3.13,  and  we  conclude  this  chapter  with  a  discussion  of  all  its  blocks,  as  we  did  for 
the  other  modules. 


Aj^,  (from  Module  D) 


Figure  3.13.  Simplified  flowchart  for  Module  C  (continued  on  the  following  page). 
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The  algorithm  starts  with  the  initialization  step  discussed  in  Section  3.4.2  (Cl),  and 
finds  the  initial  values  for  and  vj(4>m)  (C2,  C3).  Then  the  distance  )  between 

the  MaRV  and  interceptor  loci  is  found  in  Block  C4.  If  that  distance  exceeds  the  acceptance 
threshold  e  (C5),  a  correction  AA tj  is  computed  in  Block  C6,  all  relevant  variables  are 
updated  in  C7  and  C8,  and  a  new  distance  A*  is  computed  in  C3  and  C4.  Once  the  loci 
have  been  joined  at  Q*,  the  banking  time  t*BM  taken  by  M  to  reach  Q*  from  Xm(^m) 
(starting  position  for  the  turn)  is  calculated  in  C9.  Then  the  angular  distance  4>i{^*bm) 
travelled  by  I  in  time  t*BM  is  found  in  CIO  and  the  corresponding  value  for  the  turn  radius 
ri{4>i)  is  found  in  Cll. 

If  )  is  ^ss  than  the  angular  displacement  required  by  M  to  reach  Q*,  the  miss 

distance  is  zero  (C12). 

If  not,  the  miss  distance  equals 

MD  =  \Q*  -XM{m  , 

as  discussed  in  the  text  (see  Equation  (3.6)). 

Once  the  miss  distance  is  calculated,  the  result  is  fed  to  the  pK  module  in  the  battle 
manager,  as  discussed  in  Chapter  1. 

3.4.5.  MAXIMIZING  THE  PROBABILITY  OF  KILL 

Recall  from  Equation  (3.1)  that  pk  depends  not  only  on  the  miss  distance,  but  also 
on  the  lethality  radius  of  the  interceptor.  For  some  warheads,  lethality  itself  depends 
upon  atmospheric  density,  hence  upon  the  altitude  where  the  weapon  is  used.  Minimizing 
only  the  miss  distance  will  not  necessarily  yield  a  maximum  pk  for  such  weapons,  and 
the  lethality  function  or  kill  radius  must  be  inserted  into  the  optimization  loop,  as  shown 
in  Figure  3.14,  which  is  an  extension  of  Figure  3.2.  Accordingly,  when  considering  such 
syst  ems,  the  flowchart  of  Figure  3.13  is  modified  to  output  pk  values  instead  of  miss  distance 
values. 


Figure  3.14.  For  interceptors  whose  lethality  depend  upon  height  (h),  the  MD  algo¬ 
rithm  includes  the  lethal  radius  (R(h))  in  the  optimization  loop. 
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Chapter  4 

EXPERIMENTAL  RESULTS 


The  MD  Algorithm  was  programmed  in  VAX  Fortran  to  allow  comprehensive  paramet¬ 
ric  investigations,  and  in  this  chapter  we  report  the  results  of  various  computational  runs. 
The  input-output  diagram  of  Figure  4.1  illustrates  all  the  user-adjustable  parameters  of 
MD.  Of  these,  only  a  few  key  parameters  are  explored  in  this  chapter  (refer  to  the  text  to 
review  how  these  are  formally  related). 


0  = 


H0:  Maneuver/atmospheric  boundary 
hk0:  Keep-out  zone  altitude 

Interceptor  thrust  parameter  for  type  a. 
a:  Two  interceptor  types,  a  =  1  and  a  =  2 
kma:  Interceptor  drag  constants 
kDM:  MaRV  drag  constant 
mM,  mja:  MaRV  and  interceptor  masses 
Lift  coefficients 
Sj:  Effective  vehicle  area 
p  :  Atmospheric  density 
8;  Interceptor  system  time  delay 
C0V(XM>  ®):  Covariance  of  MaRV  state  and  parameter  vector 
R(y,  6)!  Orientation  of  MaRV  maneuvering  plane 
X?ia  s  Interceptor  location 

Various  numerical  parameters 


Figure  4.1.  Several  model  parameters  are  accessible  by  the  user. 
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The  parameters  discussed  in  the  various  sections  of  this  chapter  are: 


4.1. 

H0: 

4.2. 

hko' 

4.3. 

a: 

4.4. 

mI«: 

4.5. 

6: 

4.6. 

X°j: 

4.7. 

Atmospheric  or  maneuver  boundary. 
Keep-out  zone  altitude. 

Interceptor  type. 

Interceptor  mass. 

Interceptor  system  time  delay. 
Interceptor  location. 

MaRV  speed. 


Although  the  experimental  runs  were  made  by  varying  only  one  parameter  at  a  time, 
significant  effort  was  made  to  fix  the  input  vector  X m  and  the  remaining  parameters 
at  those  values  most  likely  to  reveal  important  sensitivities  and  parameter  dependencies. 
Additional  efforts  were  made  in  some  cases,  particularly  in  Section  4.1  below,  to  choose 
parameter  values  for  which  three  opportunity  windows  exist.  This  was  done  essentially  to 
reduce  the  number  of  computer  runs,  and  has  no  particular  physical  significance.  As  we 
shall  see,  a  typical  setting  for  X ^  was  (-728,  287,  10),  expressed  as  a  vector  in  units  of 
thousands  of  feet  (kft).  The  MaRV  aimpoint,  determined  by  and  the  initial  velocity 
vector  was  (125,  0,  0)  kft,  and  the  initial  speed  \V^\  was  22  kft/sec. 

The  final  speed  \Vj^\  was  2630  ft/sec,  nominal  MaRV  mass  was  15.25  slugs,  the  keep- 
out  zone  altitude  was  5000  ft,  the  nominal  interceptor  (Type  1)  mass  was  33.59  slugs,  the 
maneuvering  limit  (Ho)  was  70  kft  and  the  system  time  delay  (6)  was  0.2  sec.  No  statistical 
errors  were  propagated  in  our  tests,  and  the  covariance  matrix  cov(Xm,  9)  was  not  used 
in  our  deterministic  analysis. 

To  properly  interpret  all  the  plots  in  this  chapter,  recall  that  an  “opportunity  window” 
is  a  subset  of  waiting  times  tw  at  which  the  miss  distance  achieves  a  minimum.  Since  our 
miss  distance  predictions  are  rooted  in  the  most  recent  update  of  the  MaRV  state  Xm,  such 
windows  are  computed  for  each  such  update  at  a  rate  consistent,  with  the  battle  manager’s 
“real-time”  clock,  whose  increment  is  typically  0.1  sec.  Observe  also  that  the  occurrence  of 
triple  windows  is  rare.  Such  a  situation  typically  implies  that  the  MaRV  altitude  is  large 
(in  fact  exceeds  Hq),  and  that  it  is  currently  aimed  overhead  from  the  interceptor’s  point 
of  view.  We  shall  see  that,  as  the  MaRV  approaches  its  target,  Window  1  disappears  as 
X m  enters  the  maneuvering  region,  and  Window  3  disappears  because  tailchasing  becomes 
increasingly  difficult,  unless  the  separation  between  the  interceptor  base  and  the  MaRV’s 
target  is  large. 

Comparing  the  results  in  this  chapter  with  the  work  of  others  [13,14-18]  we  find  gen- 
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era]  agreement  in  the  conclusions,  but  have  also  discovered  some  omissions  in  our  work, 
particularly  with  reference  to  the  “field  of  view”  of  both  vehicles.  When  interceptors  are 
equipped  with  internal  homing  devices,  Or  MaRVs  with  internal  missile  avoidance  devices, 
limitations  in  the  field  of  view  can  cause  severe  errors  in  their  range  measurements.  Such 
errors  can  severely  distort  miss  distance  predictions  [14,15]. 


4.1.  Effects  of  Varying  the  Atmospheric/Maneuvering  Boundary  Ho 

Figure  4.2  illustrates  a  particular  situation  for  which  the  current  state  estimate  Xm 
produces  three  windows.  Window  i  extends  roughly  from  0  to  4  sec  and  indicates  that,  if 
an  interceptor  is  launched  now  ( iw  =  0),  the  miss  distance  will  be  very  small;  essentially  the 
same  as  if  the  MaRV  were  in  fact  a  Ballistic  Re-entry  Vehicle  (BRV),  since  maneuvering  is 
practically  excluded  in  Window  1.  The  width  of  the  window  (about  4  sec)  shows  that  the 
predicted  intercept  point  is  well  above  Hq  and  that  interception  would  still  occur  above 
Hq  even  if  interceptor  launch  took  place  4  set  later.  In  any  case,  the  best  strategy  is  to 
launch  at  the  earliest  time  at  which  Window  1  opens. 

Even  though  the  existence  of  two  additional  windows  is  mostly  of  theoretical  interest, 
since  the  policy  is  to  launch  as  soon  as  the  first  window  opens,  observe  that  a  2-sec 
tailchase  window  (Window  3)  exists.  This  usually  indicates,  based  on  the  most  recent 
state  observation,  that  the  MaRV  is  headed  Overhead  towards  a  target  significantly  far 
away  from  the  interceptor  site.  Also  note  that,  if  no  interceptor  were  fired  for  some  reason, 
the  miss  distance  would  be  about  3500  ft  if  the  interceptor  were  launched  at  the  optimal 
time  in  Region  2  ( t rs  31  sec).  s 

When  the  maneuvering  or  atmospheric  ceiling  is  lowered  from  100  kft  to  70  kft  (Fig¬ 
ure  4.3),  Window  1  is  significantly  wider,  as  would  be  expected,  since  the  MaRV  behaves 
ballistically  for  a  considerably  larger  period  of  time.  Observe  also  that  the  tailchase  window 
is  much  narrower,  down  to  about  0.25  sec.  This  is  because  the  MaRV  starts  decelerating 
much  later  in  its  trajectory  thereby  maintaining  its  initial  speed  much  longer,  and  thus 
winning  the  tailchase  during  the  first  2  sec  of  the  tailchase  region  (from  tw  —  29  sec  to 
tw  -  31  sec). 

4.2.  Varying  the  Keep-Out  Zone 

The  keep-out  zone  can  significantly  affect  the  minimum  miss  distance  achievable, 
for  two  major  reasons.  The  first  is  associated  with  the  type  of  interceptor  used,  and 
will  be  discussed  in  the  next  section.  The  second  arises  in  the  tailchase  region  (Win¬ 
dow  3).  Whenever  Window  3  is  open,  reducing  the  keep-out  height  h j.0  provides  increased 
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maneuver 


Figure  4.3.  Opportunity  windows  when  Hq  =  70  kft  (twait  in  sec  and  distances  in  104  ft). 

=  (—728k,  287k,  10k),  aim  point  =  (125k,  0,  0),  V ^  =  22  kfps,  VM^naj  = 
2630  fps,  MaRV  mass  =  15.25  slugs,  X®  =  (—20k,  0,  0),  Type  1,  int.  mass  = 
33.59  slugs,  hk0  =  5  kft,  8  =  0.2  s. 


if  h-ko  is  opportunity  for  the  interceptor  since  tailchasing  typically  arises  at  low  altitudes 
where  the  MaRV  speed  is  relatively  low  and  decreasing,  and  the  interceptor’s  speed  is  still 
increasing.  Even  reduced  by  only  a  few  thousand  feet,  the  effect  is  quite  noticeable,  as 
evidenced  by  Figures  4.4  and  4.5,  where  a  window  width  increase  of  about  1  sec  may  be 
observed.  We  emphasize  that  Window  3  closes,  not  because  the  MaRV  wins  the  tailchase 
once  again,  but  only  because  the  engagement — more  exactly  the  predicted  intercept  point 
Xq — occurs  inside  the  keep-out  zone. 

The  effects  of  varying  hko  also  show  clearly  in  Region  2,  as  shown  in  Figures  4.6  and 
4.7.  When  h *.0  is  5000  ft,  the  minimum  miss  distance  MD(t*w)  is  about  1000  ft,  whereas 
MD{t*w)  is  only  450  ft  when  hka  is  lowered  to  1000  ft.  This  is  mostly  due  to  the  fact  that 
the  MaRV  speed  is  lower  at  intercept,  and  this  is  important  considering  the  interceptor 
time  delay  of  0.2  sec  used  in  this  example. 

4.3.  Different  Interceptor  Types 

Recall  from  Chapter  2  that  a  Type  1  interceptor  is  designed  to  have  a  long  reach  and 
is  intended  for  “overlay  interception.”  In  contrast,  a  Type  2  interceptor  typically  produces 
higher  thrust,  but  for  a  shorter  period  of  time,  thereby  reaching  a  given  speed  considerably 
earlier  than  a  Type  1.  It  may  thus  be  considered  a  “last  ditch  defense”  vehicle  intended 
for  “underlay  interception”  near  the  keep-out  zone. 

Our  two  interceptor  types  are  sufficiently  different  in  that  respect  that  it  was  simple 
to  construct  a  scenario  where,  due  to  the  keep-out  zone  constraint,  a  Type  1  vehicle  was 
unable  to  intercept  a  MaRV,  whereas  a  Type  2  had  a  Window  2  extending  over  several 
sec. 


Referring  to  Figure  4.8,  the  miss  distance  remains  at  10,000  ft  for  the  Type  1,  whereas 
it  drops  to  a  minimum  of  1200  ft  in  Region  2  for  a  Type  2. 

4.4.  Varying  Interceptor  Mass 

Because  the  minimum  radius  of  curvature  of  a  maneuver  is  directly  proportional  to 
the  mass  of  a  vehicle  (see  Chapter  2),  and  the  path  length  along  the  maneuver  is  also 
proportional  to  that  radius  (see  Chapter  3),  one  would  expect  that  the  miss  distance 
would  also  be  roughly  proportional  to  interceptor  mass.  This  is  in  fact  typically  the  case, 
as  shown  in  Figure  4.9,  where  we  varied  the  interceptor  mass  from  3.36  slugs  to  335.9  slugs. 
A  converse  relationship  obviously  holds  for  the  MaRV.  As  its  weight  is  increased,  the  miss 
distance  decreases  since  its  maneuvering  capability  decreases  proportionately,  but  this  test 
was  not  executed  on  the  computer. 
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keepout  =  5000  ft. 


Figure  4.4.  The  three  windows  when  the  keep-out  height  hk0  is  5000  ft  (fwait  in  sec 
and  distances  in  103  ft). 

=  (-728k,  287k,  10k),  aim  point  =  (125k,  0,  0),  V ^  =  22  kfps, 
^Mfinal  =  ^630  fps>  MaRV  mass  =  15.25  slugs,  Xj  =  (—20k,  0,  0), 
Type  1,  int.  mass  =  33.59  slugs,  Hq  =  100  kft,  6  =  0.2  s. 


keepout  =  1000  ft 


0000 1  *  isipssiuj 


The  three  windows  when  the  keep-out  height  is  hko  is  1000  ft  (twait  in 
sec  and  distances  in  104  ft). 

X°M  =  (-728k,  287k,  10k),  aim  point  =  (125k,  0,  0),  =  22  kfps, 

VMjjnai  =  2630  fps,  MaRV  mass  =  15.25  slugs,  Xj  =  (—20k,  0,  0), 
Type  1,  int.  mass  =  33.59  slugs,  Hq  —  100  kft,  S  =  0.2  s. 


keepout  =  5000  ft. 
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Figure  4.6.  Window  2  when  h^o  is  5000  ft  (twajt  in  sec  and  distances  in  104  ft). 

XM  =  (-728k,  287k,  10k),  aim  point  =  (-10k,  0,  0),  Vft  =  22  kfps, 
Vm^ nal  =  fps,  MaRV  mass  =  15.25  slugs,  Xj  =  (5k,  0,  0),  Type  1, 

int.  mass  =  33.59  slugs,  Hq  —  70  kft,  8  =  0.2  s. 
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Figure  4.9.  The  minimum  miss  distance  in  Region  2  is  very  sensitive  to  interceptor 
terminal  mass  (Zwait  in  sec  and  distances  in  103  ft). 

X%  =  (-400k,  150k,  -50k),  aim  point  =  (0,  0,  0),  Vjh  —  22  kfps, 
^Affinal  =  5000  fps,  MaRV  mass  =  15.25  slugs,  hk0  =  5000  ft,  Xj  —  (0, 
0,  10k),  type  1  interceptor,  H0  =  70  kft,  6  =  0.2  s. 


4.5.  Accounting  for  Interceptor  System  Time  Delay 

Recently,  the  MD  Algorithm  has  been  augmented  with  the  capability  to  account  for 
interceptor  time  delays,  as  described  in  Sections  3.4.1  and  3.4.2.  These  delays  can  severely 
reduce  interceptor  performance,  particularly  if  they  exceed  0.2  sec.  Figure  4.10,  for  in¬ 
stance,  illustrates  the  miss  distance  in  Region  2.  As  the  time  delay  8  is  increased  from  0  to 
0.5  sec,  the  minimum  miss  distance  varies  from  0  to  about  4000  ft;  a  significant  spread.  For 
practical  systems,  8  ranges  between  0.1  and  0.2  sec,  and  the  corresponding  miss  distance 
ranges  between  750  and  1500  ft. 

To  show  the  influence  on  the  kill  probability  P*.,  a  lethality  table  was  folded  in,  and 
the  resulting  Pj.  values  were  plotted  in  Figure  4.11.  For  this  particular  payload  lethality 
assumption,  a  delay  of  up  to  0.2  sec  is  quite  tolerable,  but  performance  drags  off  sharply 
as  8  is  increased  beyond  that  value. 


4.6.  Interceptor  Location 

Interceptor  base  location  strongly  influences  interceptor  performance  only  when  the 
interceptor  is  launched  far  from  the  MaRV  aimpoint  or  destination.  In  two  extreme 
cases  shown  in  Figure  4.12,  where  the  most  recent  MaRV  position  estimate  was  (—728  kft, 
287  kft,  10  kft)  and  the  interceptor  was  located  at  (—400  kft,  0,  0)  and  (+400  kft,  0,  0),  the 
interceptor  missed  the  MaRV  because  it  was  unreachable  in  both  cases.  In  the  first  case, 
and  even  though  the  interceptor  was  located  400  kft  towards  the  MaRV  from  the  origin, 
the  MaRV  aimpoint  was  so  far  beyond  the  origin  that  it  was  totally  unreachable,  passing 
overhead.  In  the  second  case,  every  point  along  the  MaRV  trajectory  was  unreachable. 

When  the  interceptor  was  moved  to  (  —  100  kft,  0,  0),  the  interceptor  acquired  a  large 
exoatmospheric  window  of  about  20  sec  during  which  interception  was  assured. 

When  the  interceptor  was  moved  to  the  origin,  additional  windows  opened  up,  and 
even  a  tailchase  opportunity  showed  up,  as  seen  from  the  top  right  diagram  in  Figure  4.12. 

When  the  interceptor  was  located  at  (  —  100  kft,  0,  0.),  a  wide  exoatmospheric  window 
appeared,  as  shown  in  the  bottom  left  diagram.  When  located  at  (+100  kft,  0,  0),  there 
was  no  Window  1,  but  a  good  opportunity  in  Region  2  appeared,  together  with  a  narrow 
tailchase  window  where  the  interceptor  had  the  advantage.  This  window  was  essentially 
due  to  the  fact  that  the  MaRV  was  predicted  to  pass  overhead,  towards  its  extreme  down 
range  destination  of  (125  kft,  0,  0). 
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delta 


sensitive  to  interceptor  time  delay: 


4.7.  Varying  the  Initial  MaRV  Speed 

The  last  parameter  varied  in  our  study  was  the  MaRV  speed  |V^|.  Recall  that  this  is 
the  most  recent  speed  estimate  handed  over  by  sensors  or  radars. 

As  seen  from  Figure  4.13,  the  miss  distance  is  not  as  sensitive  to  MaRV  speed  changes 
as  it  is  to  mass  changes,  and  the  same  comment  applies  to  interceptor  speeds.  Conse¬ 
quently,  one  might  be  tempted  to  trade-off  speed  for  mass,  but  be  aware  that  the  reduction 
in  payload  lethality  associated  with  a  reduction  in  payload  mass  could  entirely  offset  any 
maneuvering  advantage  acquired.  In  some  cases  a  reduction  in  payload  mass  may  in  fact 
reduce  interceptor  performance,  but  this  possibility  is  strongly  warhead-dependent. 
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vm0  =  k  *vm00:  k  =  0.6,  0.8,  1.0,  1.2,  1.4;  vm00  =  22000 
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Figure  4.13.  The  miss  distance  (in  102  ft)  is  not  very  sensitive  to  initial  MaRV  speed 
(in  104  ft /sec). 

X*h  =  (—400k,  150k,  —50k),  aim  point  =  (0,  0,  0),  =  5000 

fps,  MaRV  mass  =  15.25  slugs,  hka  —  5000  ft,  int.  mass  =  33.59  slugs, 
Xj  =  (0,  0,  10k),  type  1  interceptor,  Hq  =  70  kft  8  =  0.2  s. 


Chapter  5 

CONCLUSIONS  AND  FUTURE  WORK 

In  allocating  interceptors  to  decoyed  MaRVs  and  BRVs,  the  battle  manager  must  care¬ 
fully  process  the  discrimination  information  handed  down  by  various  sensors,  and  must  be 
able  to  predict — or  at.  least  to  estimate — the  outcome  of  an  engagement,  sometimes  far  in 
advance  of  launch  time.  The  battle  manager  discussed  in  this  report  satisfies  such  require¬ 
ments.  Discrimination  information  is  processed  in  the  form  of  a  conditional  probability 
that  the  threat  considered  is  in  fact  a  MaRV  or  a  BRV,  conditional  to  the  observations, 
and  this  information  itself  conditions  all  the  major  decisions  made  in  the  battle  manager. 

The  outcome  of  an  engagement  is  judged  by  the  expected  asset  loss  incurred  as  a  result 
of  MaRV  and  BRV  leakage.  Because  the  interactions  behveen  individual  engagements  are 
considerably  weaker  than  interactions  within  individual  engagements,  we  have  emphasized 
one-on-one  engagements  in  this  report.  Their  outcome  is  judged  by  the  probability  of 
kill  achievable  by  an  interceptor  allocated  to  a  given  MaRV  (or  BRV).  This  probability 
is  derived  from  a  statement  about  the  lethality  of  the  interceptor  and  the  miss  distance 
achievable  by  the  interceptor.  We  have  intentionally  left  lethality  discussions  for  future 
reports,  and  have  concentrated  on  developing  a  method  for  predicting  miss  distances. 

The  algorithm  we  developed  for  predicting  miss  distances — called  the  MD  Algorithm — 
uses  geometric  and  kinematic  information  about  the  interceptor  and  the  MaRV  to  produce 
a  “best”  estimate  of  miss  distance  from  the  most  recent  MaRV  state  update  and  from 
technological  parameters  and  constraints.  In  predicting  miss  distances,  MD  accounts  for 
information  and  response  time  delays  between  threat  objects,  sensors,  and  interceptors.  It 
also  contains  a  sensitivity  framework  with  which  the  propagation  of  errors  and  uncertainties 
from  data  inputs  to  output  decisions  can  be  evaluated,  but  this  capability  has  not  yet  been 
implemented  in  the  software. 

The  MD  algorithm  was  implemented  and  thoroughly  tested,  and  several  conclusions 
became  apparent.  For  typical  system  time  delays  of  0.1  sec  and  for  nominal  MaRV  and 
interceptor  capabilities  (see  Chapter  4),  the  minimum  miss  distance  achievable  for  “late” 
engagements  (those  occurring  below  80,000  feet  altitude)  was  approximately  600  feet.  For 
“early”  engagements  MaRV  maneuvering  was  limited  and,  if  the  MaRV  was  reachable  by 
the  interceptor,  miss  distances  were  less  than  20  feet. 

Another  observation  relates  to  the  relative  sensitivity  of  miss  distance  to  interceptor 
velocity  and  maneuvering  capabilities.  It  was  found  that  the  miss  distance  was  considerably 
more  sensitive  to  the  turning  capability  of  the  vehicle  than  to  vehicle  speed.  In  fact,  miss 
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distance  was  relatively  insensitive  to  MaRV  speed.  This  suggests  that  current  development 
efforts  to  improve  vehicle  maneuverability  could  be  more  effective  than  those  to  improve 
their  speed. 

Several  things  remain  to  be  done.  Foremost  among  these  is  the  thorough  testing  of  the 
battle  manager  as  a  whole.  While  its  most  important  module — the  MD  Algorithm — was 
carefully  tested  as  an  individual  unit,  the  battle  manager  has  not  yet  been  completely 
tested  with  that  unit  in  place.  The  MD  algorithm  involves  some  heuristics,  and  several 
parameters  must  be  set  while  it  is  operating  as  an  integral  part  of  the  battle  manager. 

Such  testing  may  also  suggest  ways  to  improve  any  of  the  algorithms  discussed  in 
this  report.  For  instance,  the  straight-line  MaRV  trajectory  prediction  used  in  MD  might 
be  constrained  or  improved  by  using  additional  information  about  the  destination  of  the 
MaRV.  Further  simulation  may  also  suggest  that  the  maneuvering  test  should  be  improved, 
but  we  believe  that  this  is  unlikely  for  the  types  of  objects  discussed  in  this  report. 
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ments  to  our  work  were  suggested.  Finally,  we  would  like  to  thank  F.  McFarland  for  her 
text  editing  and  production  work,  D.  Gallant  for  her  graphics,  and  W.  Clements  for  his 
technical  editing  of  this  report. 
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APPENDIX  A 


Coordinate  Systems  and  Transformations 

A.l.  Global  (Reference)  Coordinate  System  ( Re ) 

The  global  coordinate  system  is  the  earth  center  frame  as  shown  in  Figure  A.l  [Al]. 
The  earth  center  frame  is  a  rectangular  Cartesian  frame  whose  origin  is  centered  at  earth 
center,  the  X-axis  (Xe)  is  determined  by  the  earth  center  and  the  point  0°  latitude  and  0° 
longitude,  the  Z-axis  ( Ze)  lies  on  the  equatorial  plane  and  at  90°  from  the  s-axis.  Finally, 
the  F-axis  ( Ye )  connects  the  earth  center  to  the  north  pole. 


Ye 


0  Long. 


Figure  A.l.  Specifications  of  global  coordinate  system  and  local  level  frame. 


A. 2.  Local  (Reference)  Coordinate  System  (12-System) 

As  seen  in  Figure  A.l,  the  x  and  z  axes  df  this  system  lie  on  a  plane  tangent  to  the 
earth  at  a  point  that  is  “arbitrarily”  defined  as  the  origin  of  the  local  coordinate  system. 
The  x-axis  in  this  12-system  always  points  south  and  aligns  with  the  longitudinal  meridian. 
The  y-axis  points  up,  and  the  2-axis  points  west. 

The  origin  is  usually  chosen  to  be  the  centroid  of  the  asset  locations. 


A. 3.  MaRV  Coordinate  System  (Af -System) 

The  Af -system  is  centered  at  the  MaRV  position  vector  Xm(^M )>  expressed  in  the 
12-system. 


If  ey  is  a  unit  vector  along  the  y-axis  (of  the  12-system),  and  ey^o  is  a  unit  vector  along 
the  MaRV  velocity  vector  then  the  unit  vector  triple  (eXM ,  e^,  eZM )  along  the  MaRV 
axis  triple  y^,  zm)  is  defined  as  follows: 


—  €y  o 
VM 


eyM^w)  —  x  e 


XM  ? 


Ym  x  eg 
|vj  X  e. 


(Fj»  x  e,)-e,  >»AVJx  AX ,(iw)  =  0 


r  {v** x  €x)'ey  - 0  A  V°M  x  AXl{tw)  =  0 

£ fnt)  x  Sr  (Ajr'(M *  v°m)  ■«.  > » 


where  e^X/  is  a  unit  vector  along  the  interceptor  launch  direction  vector  A Xj.  Note  that 
for  each  new  value  of  twait  we  get  a  new  coordinate  system  (xm,  VM ,  ZAf)(t wait)* 

The  purpose  for  defining  the  M- system  this  way  is  to  facilitate  the  analysis  and  compu¬ 
tation  of  MaRV  and  interceptor  maneuvers.  The  nominal  plane  in  which  these  maneuvers 
take  place  (the  nominal  “action  plane”)  is  the  (e^o,  e^X;  )-plane.  When  maneuvers  are 
not  in  that  plane,  they  are  transformed  into  that  plane  by  the  rotation  R  discussed  below. 
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In  the  definition  of  eZM ,  a  test  is  needed  to  avoid  conditions  where  zero  vectors  would 
be  chosen  as  unit  coordinate  vectors,  as  would  be  the  case  when  e&Xr  and  ev°  are  collinear, 

1  M 

i.e.,  evo  x  eAXj  = 

Another  test  is  needed  to  assure  that  eZM  is  pointing  in  the  right  direction,  and  this 
requires  calculating  the  dot  products  with  ey ,  as  shown. 

To  minimize  computation  time  in  determining  the  M-coordi nates  as  the  engagement 
unfolds,  we  use  the  following  flowchart  (Figure  A. 2). 


Figure  A. 2.  Determining  the  MaRV  coordinate  system. 
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A. 4.  Coordinate  Transformations 

Transformations  are  needed  in  the  miss  distance  calculations  because  vehicle  height 
(altitude)  in  the  R- system  is  an  important  atmospheric  parameter.  Altitude  is  simply  the 
projection  Vrv  of  the  position  vector  Vr  onto  the  y-coordinate  and,  because  the  curvature 
radius  of  the  MaRV  depends  upon  its  altitude,  which  is  expressed  in  the  iZ-system,  we 
need  to  transform  M- points  into  R- points  and  vice-versa.  These  transformations  require 
both  a  translation  and  a  rotation.  Consider  Figure  A. 3  where  a  transformation  from  M 
to  R  is  illustrated. 


Figure  A. 3.  Transforming  a  vector  Vm  in  the  M-c oordinate  system  to  a  vector  V  in 
the  R- system. 

Rotations  are  done  in  three  steps  [A2]. 

1*  Rotate  Vju  about  the  zj^- axis  by  the  (elevation)  angle  6  between  x m  and 
the  x-z  plane  (sin(0)  =  —  eyM  *ey),  yielding  a  rotation  matrix 

cos  6  sin  0  0  ‘ 

RZm(9)  —  —  sin  6  cos  6  0  .  (A4) 

.  0  0  1. 
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2.  Rotate  the  rotated  vector  RZm(0)’Vm  about  the  yj^  axis  by  an  angle  (3 
between  xm  and  the  x-y  plane,  yielding  a  rotation  matrix 


RvmW)  ~ 


cos/?  0  sin  f3 
0  1  0 
—  sin  f3  0  cos  (3 


(A5) 


3.  Rotate  the  rotated  vector  RyM((3)'  RZm((3)’Vm  about  the  xm  axis  by  an 
angle  7  between  zm  and  the  x-y  plane,  yielding  a  rotation  matrix 


Rzm( 7)  =  |  0  cos7  sill  7 

.  0  —  sm  7  cos  7  j 

The  composite  rotation  of  Vm  thus  yields  a  vector 


1 


0 


0 


(A6) 


Vr  =  RZM(j)'RyM(/3)'RXM(0)VM 


(A7) 


in  the  R-system. 

Since  the  origins  of  the  M- system  and  the  R-system  are  different,  we  also  need  to 
translate  Vm-  To  be  consistent  with  the  meaning  of  the  rotation  angles,  this  must  be  done 
first,  before  rotating  Vm-  The  overall  transformation  is  thus 

VR  =  T(VM)  =  RzM(l)-RyM(P)'RxM(0)lvM  -  Xm]  ,  (A8) 

where  Xm  is  the  position  vector  of  the  M- system  in  the  R-system.  Thus 

Vrv  =  Projy[RZM(j)-RyM(l3)>RXM(0){VM  -  XM]}  -  (A9) 
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APPENDIX  B 


Simplified  Equations  of  IVIotion  for  Interceptor  Type  I  (High  Altitude) 
B.l.  Mass  (Figure  B.l) 


weight 


Figure  B.l.  Weight  profile  of  a  type  I  interceptor. 


+  Fi(<)  *  +  Zi(t )  *  rnji$(t) 

Xi(<)  =  1  +  sign(1-61"<) 

2 


^1  +  sign(3.83  —  /)^  ^1  +  sign(/  —  1.61)^ 

1  +  sign (2  -  3.83)  \ 

2  J  ’ 

mIn(t)  =  230.5(1  -  0.449/)  +  ( m°n  -  33.59)  ,  0  <  /  <  1.61  , 

m/i2(*)  =  85.91(1  -  0.159/)  +  (m°n  -  33.59)  ,  1.61  </<  3.83  , 

ra/i3(/)  =  rn°n  ,  3.83  <  /  . 


(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(B7) 
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B.2.  Thrust  (Figure  B.2) 


Figure  B.2.  Thrust  profile  of  a  type  I  interceptor. 


Tn(t)  =  Xi(0  *  Tm(t)  +  Yi(<)  *  TI12(t)  +  Zj(<)  *  Tm(t)  ,  (B8) 

Tm(t)  =  625,500  ,  0  <  t  <  1.61  ,  (B9) 

Tn2(t)  =  114,000  ,  1.61  <<<3.83  ,  (BIO) 

2j13(<)  =  0  ,  3.83  <  t  .  (Bll) 


B.3.  Acceleration  Due  to  Thrust 

aIi(t )  =  ^0(0  *  aIll(t)  +  Fj(<)  *  a/i2(<)  +  Z\ (t)  *  a/i3(<)  , 


ajn(t)  = 

o/12(0  = 


625,500 


230.5(1  -  0.449t)  +  (m^  -  33.59) 

_ 114,000 _ 

85.91(1  -  0.1590  +  (m^  -  33.59) 


a/i3(0  =  0  ,  3.83  <  t  . 


0  <  t  <  1.61  , 

1.61  <  t  <  3.83 


(B12) 

(B13) 

(B14) 

(B15) 
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B.4.  Speed  [B1,B2] 

Vn(t)  =  Xi(t),Vm(t)  +  Y1(t)*VIn(i)  +  Z1(t)*Vm(t)  , 

/k  k 

- —dx  —  —  In  (a  +  bx )  , 

a  +  bx  b 


Vjn(t ):  «  =  +  196.9  ,  6  =  -103.5  ,  -  =  -6043  . 

b 


(B16) 

(B17) 

(BIS) 


V,n (1)  =  -60431n(m5i  +  196.9  -  103.5r) 

= -6043(ln(m5,  +  196.9  -  103.51)  -  lnfm?!  +  196.9)]  ,  0  <  t  <  1.61  .  (B19) 


Vn,U):  a  =  m$,  +  52.32  ,  b  =  -13.66  ,  -  =  ■  83-16  . 

b 


VJ12(t)  =  -83461n(m5j  +  52.32  -  13.66r)  +  F/n(1.61) 

r=1.61 


(B20) 


=  -6043(111(7775!  +  30.26)  -  ln^i  +  196.9)]  +  83461n(ra5i  +  30.33) 

-  83461n(7775i  +  52.32  -  13.66*)  ,  1.61  <  <  <  3.83  .  (B21) 

Vm{t)  =  23031n(7775i  +  30.30)  +  6043  ln^j  +  196.9) 

-  83461n(Tn5i  +  52.32  -  13.66*)  ,  1.61  <  *  <  3.83  .  (B22) 

.  .  f  Vln  I1  -  (tFtP)]  .  3.83  <  1  <  3.83  +  kDh 

Vm(<)=<  1  /J  ,  (B23) 

.  0  ,  3.83  +  fcp/j  <  * 

where 

1/13  =  V/i2(3.83)  =  23031n(ra5!  +  30.30)  +  60431n(7775!  +  196.9)  —  83461n(7775!) 


V?  f 

'713  [ 

t— 3.83)  \  1 

- 1 

^  kDI  1  )\ 

.0, 

(B21) 

(B22) 

(B23) 


B.5.  Position  [B1,B2] 

Xn (*)  =  *i(<)  *  Xni(t)  +  Fi(*)  *  XIn(t)  +  Zx{t)  *  XIn(t)  , 
Fact:  J  In  (a  a:  +  b)  =  ^aa:  ^  in(a;c  _|_  _  x  ? 

X/n(*):  a  =  —103.5  ,  b  =  mjj  +  196.9  . 


(B24) 

(B25) 

(B26) 
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Xm(t) 


Xn2{t):  a 

Xn  2(0 


Xxuit) 


6043  1111(771?!  +  196.9)*  +  (TO/1  +  ln(™?i  +  196-9  “  103-5r)  +  r[  J 


6043<f  (In  (twJi  +  196.9)  +  1  )  *  +  (mJi  +  19^-9  103-5<)  +  196.9  _  103.5 *) 

S  lUo.O 


/m°n  +  196.9  \  ,  ,  o  _  rt,l 
“  ( — 103.5 - J  +  196-9)} 

6043|ln(?7i?j  +  196.9)  +  1  ]*  +  +  1-90  -  *^  111(771?!  +  196.9  -  103.5*) 

-  (m/\o3^5  "  ~)  ln(mJi  +  196-9) |  ,  0  <  *  <  1.61  . 

-13.66  ,  6  =  m?  ,  +  52.32  . 

X/n(1.61)+  /  VI12(r)dr 

t 

Xjn(1.61)  +  [23031n(m?i  +  30.30)  +  60431n(m?i  +  196.9)]r 

1 .61 

/ 771?.  +  52.32  —  13.66r\  0 

(  +13  66 - )  ln^m”1  +  52'32  _  13*66r)  +  r 


+  8346 


X/ii(1.61)  +  [2303  ln(T7i?i  +  30.30)  +  60431n(m?i  +  196.9)  +  8346]* 
-  [37081n(m?i  +  30.30)  +  9729\n(m°n  +  13,437] 

+  8346  [  ( +  3.83  -  f )  ln(77i?i  +  52.32  -  13.66*) 

L  \13.66  / 

_  f  13  (}6  ^n(m/i  +  30.30)j  , 

X/ii(1.61)  +  [2303  ln(  77i?i  +  30.30)  +  60431n(ro?i  +  196.9)  +  8346]* 

+  8346  ( -^IL-  +  3.83  -  A  ln(m?i  +  52.32  -  13.66*) 

\  13.66  / 

-(611.0t7i?j  +  22, 238)  ln(m?i  +30.30) 

+  9729  ln(m?i  +  196.9)  +  13,437  ,  1.61  <*<  3.83  . 


(B27) 

(B28) 


(B29) 


(B30) 
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Xm (t)  =  X/12(3.83)  +  f  VI13 (r)dr  =  X7l2(3.83)  +  Vfn  \(l  +  r 

J  r=3.83  LV  kDI  1/ 

=  *»<«*)  +  V?13  [(l  +  H)  t  -  ^  -  3.83  (l  +  |5L)]  , 


3.83  <  t  <  3.83  +  kDI1 


2kD11 


3.83 


Xm(i)  =  X/i3(3.83  +  k^n)  ,  3.83  +  <  t 


(B31) 
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APPENDIX  C 


Simplified  Equations  of  Motion  for  Interceptor  Type  II  (Low  Altitude) 
C.l.  Weight/Mass  (Figure  C.l) 


weight 


0  ti  =  1  t2  =  2  time  t  (sec) 


Figure  C.l.  Weight  profile  of  a  short-range  type  II  interceptor. 

mj2(t)  *  mIn{t)  +  y2(<)  *  mj22(t)  +  Z^  *  mm{t)  , 

(Cl) 

X2(t)  = 

1  +  sign(  1  -  t) 

2 

(C2) 

n(<)  =  | 

f  1  +  sign  (2  ~t)\  /I  +  sign  (<  -  1)\ 

(  2  A  2  )  ' 

(C3) 

W)  = 

^1  +  sign  (t  —  2)^ 

(04) 

mm(t)  =218.8(1  - 

-  0.714<)  +  m°I2  -  15.62  ,  0  <  t  <  1  , 

(05) 

mj22(t)  =109.4(1  - 

-  0.429<)  +  mj 2  -  15.62  ,  1  <  t  <  2  , 

(06) 

mm{t)  =m°j 2  , 

2  <  t  . 

(C7) 
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C.2.  Thrust  (Figure  C.2) 


Figure  C.2.  Thrust  profile  of  a  short-range  type  II  interceptor. 


Tn{t)  =  X2(t)  *  Tm{t)  +  y2(t)  *  TI22(t)  +  Z2(t)  *  Tm(t)  , 
Tm(t )  =  2  x  106  ,  0  <  t  <  1  , 

TI22{t)  =  105  ,  0  <  t  <  2  , 

7/23(0  —  0  ,  2  <  t  <  oo 


(CIO) 

(Cll) 


C.3.  Acceleration  Due  to  Thrust 


aI2(t)  —  A  2 (0  *  Oj2i(t)  +  F2(t)  *  0/22(t)  +  i?2(t)  *  0/23(0  5 

2  x  106 

=  218.8(1  -  0.7141)  +  m$2  -  15.62  ’  0  -  '  -  1  ’ 


(012) 


(013) 


0/22(0  ~ 


109.4(1  -  0.4290  +  m%  -  15.62 


,  1  <  t  <  2  , 


(C14) 


0/23(0  =  0  ,  2  <  t  . 


(C15) 
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C.4.  Velocity  (Speed)  [B1,B2] 

Vn(t)  =  X2(t)  *  VI21 (t)  +  Y2(t )  *  V122(t)  +  Z2(t)  *  VI23(t )  . 

_  t  kdx  k .  .  ,  . 

Fact:  /  - —  =  —  m(a  +  bx)  , 

J  a  +  bx  b 

k 

Vm{t):  a  =  203.18  +  m°n  ,  b  =  -156.2  ,  -  =  -12,803  . 

b 

VI21(t)  =  — 12,8031n(203.18  +  m°I2  -  156.2r)  * 

=  -12, 803[ln(203.18  +  m°J2  -  156.2t)  -  ln(203.18  +  m°l2 )]  ,  0  <  t  ; 

Vi22{t):  a  =  93.78  +  m°I2  ,  b  =  -46.98  ,  7  = -2131  , 

0 

Vm{t)  =  12, 803[ln(203.18  +  m°12 )  -  ln(46.98  +  m°I2)] 

-  2131  [ln(93.78  +  m52  -46.98)  -  ln(46.98  +  m?2)]  ,  l<t<2 

(  12,803[ln(203.18  +  m°I2)  -  ln(46.98  +  m°I2 )] 

Vm(t)  =  { 

l.  —  2131[ln(m72)  —  ln(46.98  +  mj2)] 

F7°23  =  12, 803 ln( 203.18  +  m52)  -  10, 672 ln(46.98  +  m?2) -2131  ln(m52)  , 


C.5.  Position  [B1,B2] 

XI2(t)  =  X2(t)  *  xm(t)  +  Y2(t)  *  XI22(t)  +  Z2(t)  *  XI23(t)  , 

Fact:  J  k  ln( ax  +  b)dx  =  k  ^  ln(aa:  +  b)  —  x  , 

X/21(t):  a  =  -156.2  ,  b  =  203.18  +  m°12  ,  A:  =  —12, 803  . 


Xm(t) 


=  12, 803 j 


ln(203.18  +  m52)i 


+ 


^203.] 


18  +  mQl2  —  156. 2r 


+156.2 


^  In  (203.18  +  m72  —  156.2r)  +  r 


(C16) 

(C17) 

(CIS) 

1  .(C19) 
(C20) 

•  (C21) 

(C22) 

<  t  ■ 
(C23) 

(C24) 

(C25) 

(C26) 
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=  12, 803<[ln(203.18  +  m°J2)t  +  i  +  ( 1.3  +  ^  ln(203.18  +  rn52  -  156.2/) 

[  \  156.2  ) 


ln(203.18  +  m 


12)  | 


=  12, 803  j  (In  (203.18  +  m°j2)  +  l)t-  (L3  +  ln( 203.18  +  m°I2) 

+  (1-3  +  ^^~<)ln(203*18  +  Tn/2-156-20}  »  0</<l  •  (C27) 


0  <  /  <  1  . 

Xm{t):  a  =  -46.9  ,  b  =  93.78  +  m°J2  ,  k  =  -2131  ,  (C28) 


Xm(t)  =  XI  2i(l)+  f  VI22(r)dr 

J  r—\ 

=  X/2i(l)  +  r [12, 803 ln( 203.18  +  m°I2)  -  10, 672  ln(46.98  +  m°J2)} 

nn  r /46.9t-  (93.78  +  ,  n  .  1* 

-  2131  ^ ln(93.78  +  m°I2  -  46.9t)  -  rj 


(C29) 


=  X/2i(l)+  (12, 803 In  (203.18  +  m®2)  -  10, 672 In  (46.98  +  m°I2)  +  2131)/ 

-  12, 8031n(203.18  +  m°12)  -  2131  +  (8541  -  45.44 ra?2)ln(46.98  +  m°I2) 

+  (4261  +  45.44 m°I2  -  2131/)  ln(93.78  +  m°I2  -  46.9/)  ,  1  <  /  <  2  .  (C30) 


Xm(t): 


Xm(t)  =  Xm{2)  +  jT  ^  Vm (r)dr  =  Xm{2)  +  V?23  (l  +  r  - 


n  t 


2kDJ2 


J  2 


&D/2 
2  <  /  <  2  +  &£>/2 


/2  \ 
2kon  ) 


(C31) 
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APPENDIX  D 


Some  Derivatives  for  Using  the  Newton- Raphson  Method  in  Solving  f(x)  =  0. 

In  using  the  root  solving  routine  RTSAFE  [10],  a  significant  speed-up  is  obtained  if 
the  derivatives  of  f(x )  are  provided.  In  this  appendix,  we  present  these  derivatives  for  the 
initialization  problem. 

“Solve: 

F  (^wait?  At/)  =  At/)  —  Ai2(twajt,  At/)  =  0 

for  At/,  twa it  fixed.”  (Dl) 

We  need  the  derivative  (F°)'(twa ;t,  At/)  of  F°(twa ;t,  At/). 


D.l. 

F  (twaitj  At/)  =  d^//(twait,  At/)  —  AR(twajt,  At/)  ,  (D2) 

)  (^wait5  At /)  =  (twait?  At/)  —  AR  (twajt,  At/)  .  (D3) 


D.2. 


^At/(^wait)  At/)  2[rji/(twajt,  At/)r/(twaj^,  At/)]  ^  cos(#c)  +  ^(t^ait,  At/)(sin^)  , 

(D4) 

where  0 °  is  the  colfision  angle  between  the  MaRV(M)  and  the  interceptor  (/). 


Define 


Atjtf  —  twait  "h  At/ 


(D5) 


Then 


and 


Now 


rM(^wait?  Afj)  —  ^wait  d"  At])  —  5 

r/(^wait?  Aij)  ~ 


(^wait?  A^j)  — 


d^M7  ^rM  ^M7  ^r7 
^rAf  5A</  dr  /  dAtj 


(D6) 

(D7) 


D.S. 


(D8) 
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dd°MI  (  a(h  (rM{^M)\i^  ,  .  Qo 
dr. 1  c>  V  rj(A<j)  / 


rM{i-wa.iti  At j)  V^^Atj  ^wait)  —  ^’M(Atjif)  , 


(DIO) 


drMjAtM) 

d(Ati) 

9Xm{A  tji) 
d(Atj) 


2mMexP{XM{A^M)-e>  -ij 


PoCl(c*,  Ma.ch) mSm 


(Dll) 


2raji/  exp 


XM(AtM)'e } 


-} 


poC,£(a,  Mach^^M 


dXMjAtMyty 

hod(Atj) 


)  • 


(D12) 


|Vjj(A(M)|evj,  > 

f  ivii , 


0  <  A/jif  <  to 


where 


|Vm(A<m)|=<  1fmI  t1  -  iAlMDMTo)]  »  o<t0<a/m 

T0  =  T70  (Ho-X^'Cy)  , 

*M'ev 


(D13) 


(D14) 


(D15) 


sin  0®  —  evo'ey  ,  e™ 


A  XI, 


1  m  ’ 


(Die) 


k  -2F&min{go,Xgf-e,} 

(usi  -  i»sd  (iyai + ivii)  ’ 

bjj/  is  the  speed  of  M  at  last  observation,  and  V,,  is  the  final  (impact)  velocity  of  M. 
If  fjy *ej,  '  — 10— 3,  set  y^frX!  —  ]03. 


(D17) 


ri{Atj) 


2mj(Atj)  exp  |  ^  —  1^ 

PqCl{&,  Mach)/ Sj 


(D18) 
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drj 

dKTj 


2  exp 


Xi(Ati)»ey 


Po('l(<*,  Mach)/ Sj 


dmj(Atj)  rni{Aij )  / dXj(Atj)’ey 
dAtj  ho  \  dAtj 


(D19) 


m-i(Ati)  =  X(Atj)m-n(Ati)  +  Y{  A*j)mj2(Atj)  +  Z(  Af,/)m/3(  Ati)  . 

For  the  type  I  interceptor,  the  equations  for  mji(Atj)  are  provided  in  Appendix 
mm(t)  =  230.5(1  -  0.449*)  +  (m°n  -  33.59)  ,  0  <  *  <  1.61  , 

mn2(t)  =  85.91(1  -  0.159*)  +  {mQn  -  33.59)  ,  1.61  <  *  <  3.83  , 

mI13(t)  =  m/j  ,  3.83  <  *  ( *  in  seconds,  mQn  typically  33.59  slugs)  , 

{-103  ,  0  <  Atn  <  1.61 

-13.7  ,  1.61  <  Atn  <  3-83  . 

0  ,  3.83  <  Atn 

For  the  type  II  interceptor,  the  equations  for  mn(Atn)  are  in  Appendix  B. 

m/2i(t)  =  218.8(1  —  0.714*)  +  (wijj  15.62)  ,  0  ^  ^  1  , 

mI22(t)  =  109.4(1  -  0.429*)  +  ( m°I2  -  15.62)  ,  1  <  *  <  2  , 

ra/23(*)  =  mI 2  .  »  2  <  *  , 


dmnjAtn) 

dAtn 


{-156  ,  0  <  AtI2  <  1.0 

-46.9  ,  1.0  <  AtI2  <  2-0 
0  2.0  <  AtI2 

(Note:  AtM2  =  <wait  +  A*j2.) 


^d(Atjj^  ~  I vi(Atl)\evp  >  ev°  -  eAXT  • 


(D20) 

(D21) 

(D22) 

(D23) 

(D24) 


(D25) 

(D26) 

(D27) 

(D28) 


(D29) 


[AtM 

Ai2(*wait,  A ti)  =  R°MI  -  /  \VM(r)\dr  -  /  V/(r)dr  , 

Jo  Jo 


R°MI  —  d °MC  +  5 


dAi? 

dAtj 


4/C  —  l-^Af  -  Xc\  1  <$CI  -  \%C  -  -^/l  5 

=  -[|Vrjf(A<jf)l  +  |V/(A</)|]  . 


(D30) 


(D31) 

(D32) 


(D33) 
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APPENDIX  E 


Derivatives  for  Using  DBRENT  [10]  to  Minimize  the  Function  |Arj(<^M)| 

Recall  from  the  main  text  (Chapter  3),  that  we  need  to  find  the  value  of  (f>M  at 
which  A ri((j>M)  =  \r'j(4>M)\  -  is  a  minimum.  Having  the  derivatives  of  A rj((f>M) 

considerably  speeds  up  the  process  of  finding  this  minimum  (see  Appendix  G). 

E.l. 

Ar'j((f>M)=  (|r;  ((/>M)  I)'  -  (|r/(0M)  I)'  .  (El) 

E.2. 

\A(M\=  [(M4>M)?  +  (B(4,M)f\12  , 

,1 _  dr'j  _  dr'j  dA  dr’j  dB 

1  d4*M  BA  d(j>M  dB  9<f>u 


(E2) 

(E3) 
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(A Im  fixed)  . 


Ah^j(^At^j i  4* M )  ey  *  fti^)  I rM  (  A/jj/ )  |)  X 

((sin  4>m)  exM  +  (1  —  cos  4>m)  cj/m)]  » 


R 


A 


ZMVM*M 


rM(^hl)  =  ‘I’M)]  ,  • 

1  <Pm 

(  2mI(AtI)  exp  | Cy  -  l| 

|rj(Aij)|  =  min 


poCi(a,  Mach)/Sj 


10e 


(E9) 

(E10) 

(Ell) 

(E12) 


E.4. 


9A(4>m) 

d<t>M 

drM(<j>M) 

d<f>M 


(sin  4>m)—  ^rj~~ ^  +  \tm(<I>m)\ (cos  4>m)  » 

0<PM 

\tm{4>m)\  dAflM{<f>M ) 


fio 


d(f>M 


dAh 


(E13) 

(E14) 


M 


d(j>M 


ey ' 


( T 5  ft  ■>  ^)|^m(  A/j|f) |  ((cos  4>m)  eijif  "fi  (sin  4>m)  eyjif  )]  •  (E15) 


E.5. 


d<f>M 

Note:  tm{4>m)  =  ^Af-) 


_  ^CQS  _  |rM(^,M)|  (sin^ji/) 


(E16) 


E.6. 


ri(Ati)  =  ri(Ati,  <f>j ) 


<^J=0 


(E17) 


f  2m,(Ai;)  eXp(Xl(^H£  +  At;(Ay,(fa))  _t| 

|r,(^)l  =  j - .oCiK  M«h),s,  - ’  10‘ 


(E18) 
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Ahj(Atj,  <f>i{<f>M )) 

=  \ri\  (1  -  cos  4>j  (<J>m))  £zm  x  ev£>  +  \rl\  (sin  <f>i  (<f>M))  eyf 


B(4>m)  cos  0 °  -  ^(^m)  sin  0° 
\rj(M\ 


4>i(4>m)  =  cos  1 

9\rij^M) |  _  C?A dfa{<hri 

0(4>m)  h0  d<f>i  d<f>M 

c)h 

1  -  \rl\  (is'n<f>i)eZM  x  €yo  +  (cos 4>i)evp )  ^  , 

1  dU(<f>M) 


dfa 

d<t>i<t>M )  d 


'ey  ,  (E19) 
(E20) 

(E21) 

(E22) 


d<f>M  d<f>M 


(arccos  W  (<£m)) 


U(4>m)  = 


(l  -  W  (</>m)2) 

0  <  arccos Zy  ^ 

B(M  COS  8°  -  A(<f>M)  sin  0 °  E(<£M ) 


2\1/2  d<f>M 


\t\Um)  :  ’ 

w(M  =  KlWjf)  =  [A2(«  +  b2(«W)11/2  , 

o  l  svjM  v(M_ew(M 

d<f>M  d<l>M  W(4>m)  d<f>M  W(4>m?  9(f)  m 

V(<I>m)  =  cos 0®  -  A(<f>M)  sin , 

?YMm1  =  (cos  (jj  Msk)  _  (sin  ^0)  MW  , 

d<t>M  '  d<f>M  v  '  d(j>M 


dW(M 

d<f>M 


<Hm 


(A2  (4>m)  +  52  (<^m))1/2 


(E23) 

(E24) 

(E25) 

(E26) 

(E27) 

(E28) 

(E29) 
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APPENDIX  F 


Some  Derivatives  Needed  to  Compute  the  Intercept  Point  Using  RTSAFE  [10] 


Problem 


“Solve 


F(tw ait?  tc)  —  0  tc  i  ^wait  fixed. 

F(<wait,  tc)  =  t  |F/(r)| dr  -  [ Ux(tcw  f  +  Uy(icw)2  +  Uz(tcw )2]1/2  , 
Jo 

„„  ,  f,  ,  n,\VM(r)\ir 

Uxylcw)  Mix  '  |  ’ 


Uy(tcw)  ~  X0MIy  +  / 

Jo 


r<- 

)  Wl 

r,t”  VJUFj,(t)|* 


mt  \-x°  +  f  YMMlMl 

Uzytcw)  —  MIz  *  J ^  \V]tf\  5 

i  IVmI  ’  '  o<r<-n> 

IVmWI  =  |  Wl  [> - T5?]  • 

’  r°<0:ST 


dF(*wa u,  tc) 
dtc 


_  dF  9Uj{ty^ ait;  ^c) 

c)Li1  dtc 


n  —  4 


^l(^wait5  *c) 


[  \Vi(t)\(It  , 

t/0 


^(<wait,  <c) 


5ZYi(iwait}  ^c) 


=  IWc)|  • 


W2(<c)=^(^)  , 


(F10) 
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dF 

dU2 


-A~ll2Ux{tcw) 


A  —  Ux(tcw)2  +  Uy(tcw)~  +  Uz(tcw )2 


( F 1 1 ) 


dU2 

dtc 


\VM(tcu,)\ 


(F12) 


F.4. 


^3  (4l»)  —  Uy{tcw)  J 

^  =  -A-WU,(icw) 

m_Ym_ ,,  >, 


(F13) 

(P14) 

(F15) 


F.5. 


—  Uz(tcw)  , 


(F16) 


dF 

dUA 


—A~1/2Uz(tcw) 


dU 4 
dtcw 


Ym. 

n\ 


\VM(tcw)\ 


(F17) 


F.6. 

dF(^t-wa it,  ^c) 
dtc 

A-1/2 

=  |F/(^C)|  | pro  |  +  Uy(tcw)V^y  +  Uz(tcw)V^z]  \VM{tCw)\  ,  (F18) 


-4  —  [f^x(^CTO)  +  Uy(tcw )2  +  f7z(fcu»)^]  5 


(F19) 


0  ^  ^  T0  ^ 

0  ^  Tq  ^  ^civ  ? 


r0  <  0  <  * 

cw 


(F20) 
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7 


C-e,  ’ 

|VJ&-% 

1  >  10-3  , 

) 

T0  =  < 

1  io4 , 

VI 

o 

'ev\ 

<  10-3  , 

Ho  <  X^*ey  , 

(F21) 

TT 

o 
1— 1 

1 

VI 

o 

K-e,\ 

<  1CT3  , 

Ho  >  X^*ey  , 

&DM  =  < 


{  ~2 Vm  irnnjHo,  Xlcj^} 

evM°’ey(vM  ~  Vm)(Vm  +  Vm) 
103  , 


vm^vm~^  and  ev°’ey  <  _1°  3 
otherwise 


(F22) 


F.7. 

F.7.1. 


F.7.2. 


F.7.3. 


/ 


ptCw 

/  I  vM( 

Jo 


—  Vjlftcw  —  X Mli^cw)  *  0  —  tew  ^  T0 


r)\dr  =  Vj^tq  +  Vj(l-(r-Tto)fc)dr 

Jtq 


Vmt o  +  Vmt 


tcw  VfijkT2 

tcw 

+  Vm^ot 

o 

T0  ^ 

TO 

to 


T/0  yO  jl  2 

T/0  4  .  T/0  i  4  t/0  j  2  vMKlcw  , 


—  I/° 

- 


^Tq  kt2v 
tew  +  kr()tcw  -  — - 


ctu )  ? 


0  ^  Tq  ^  tcw. 


ftew 

/  |Vj,(r)|rf 

Jo 


T  = 


AT 


kf2 

4  __  Z_CW 

Inn  ~ 


A 


(F23) 


(F24) 

(F25) 

(F26) 


(F27) 


X Mzitew)  i  To  <  0  <  tcw  ,  (F28) 


where  k  =  Jfe^ . 


Ill 


F.8.  Concluding 


9.F'(^wait,  tc 
dtc 


=  |V>((.)| 


-  A(tcw )  ^2\VM(tcw 


IX 


o 

Mix 


V° 

VM  x 


APPENDIX  G 


Newton-Raphson  Method  Using  Derivatives 
Numerical  Recipes 

Algorithm  MD  uses  several  standard  algorithms  to  find  the  root  of  a  nonlinear  function 
f(x),  i.e.,  to  solve  f(x)  =  0,  and  to  find  the  minimum  of  such  a  function.  In  this  Appendix, 
we  include  discussions  and  listings  of  the  “recipes”  we  have  implemented  in  MD.  These 
were  extracted  with  minor  modifications  from  Reference  [10]. 

G.l.  Computing  the  Root  of  f(x) 

Since  the  derivatives  of  F°(twait,  A tj)  were  derived  analytically  (see  Appendix  D), 
we  employ  the  Newton-Raphson  method  using  first  derivatives  (Function  RTSAFE)  to 
compute  its  roots. 

Newton-Raphson  Method  Using  Derivatives 

Perhaps  the  most  celebrated  of  all  one-dimensional  root-finding  routines  is  Newton’s 
method ,  also  called  the  Newton-Raphson  method.  What  distinguishes  this  method  from 
others  is  the  the  fact  that  it  requires  the  evaluation  of  both  the  function  /(«),  and  the 
derivative  f(x ),  at  arbitrary  points  x.  Geometrically,  the  Newton-Raphson  formula  con¬ 
sists  of  extending  the  tangent  line  at  a  current  point  xi  until  it  crosses  zero,  then  setting 
the  next  guess  to  the  abscissa  of  the  zero-crossing  (see  Figure  G.l).  Algebraically,  the 
method  derives  from  the  familiar  Taylor  series  expansion  of  a  function  in  the  neighborhood 
of  a  point. 


Figure  G.l. 


Newton’s  method  extrapolates  the  local  derivative  to  find  the  next  esti¬ 
mate  of  the  root.  In  this  example  it  works  well  and  converges  quadrat! cally. 
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f(x  +  8)*f(x)  +  f'(x)8  +  ^-62  +  ...  .  (Gl) 

For  small  enough  values  of  6,  and  for  well-behaved  functions,  the  terms  beyond  linear  are 
unimportant,  hence  /( x  +  8)  —  0  implies 


6  = 


/(*) 

/'(*) 


(G2) 


Newton-Raphson  is  not  restricted  to  one  dimension.  The  method  readily  generalizes  to 
multiple  dimensions,  as  we  shall  see  in  the  discussion  to  follow. 


Far  from  a  root,  where  the  higher  order  terms  in  the  series  are  important,  the  Newton- 
Raphson  formula  can  give  grossly  inaccurate,  meaningless  corrections.  For  instance,  the 
initial  guess  for  the  root  might  be  so  far  from  the  true  root  as  to  let  the  search  interval 
include  a  local  maximum  or  minimum  of  the  function.  This  can  be  death  to  the  method 
(see  Figure  G.2).  If  an  iteration  places  a  trial  guess  near  such  a  local  extremum,  so  that 
the  first  derivative  nearly  vanishes,  then  Newton-Raphson  sends  its  solution  off  to  limbo, 
with  vanishingly  small  hope  of  recovery.  Like  most  powerful  tools,  Newton-Raphson  can  be 
destructive  when  used  in  inappropriate  circumstances.  Figure  G.3  demonstrates  another 
possible  pathology. 


Figure  G.2.  Unfortunate  case  where  Newton’s  method  encounters  a  local  extremum 
and  shoots  off  to  outer  space.  Here  bracketing  bounds,  as  in  RTSAFE, 
would  save  the  day. 

Why  do  we  call  Newton-Raphson  powerful?  The  answer  lies  in  its  rate  of  convergence. 
Within  a  small  distance  e  of  x  the  function  and  its  derivative  are  approximately: 


f{x  +  e)  =  f(x)  + ef'(x) +  +  ■■■  , 

f’(x  +  t)  =  /'<*)  +  ef"(x)  +  ■■■  .  (G3) 
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Figure  G.3.  Unfortunate  case  where  Newton’s  method  enters  a  nonconvergent  cycle. 

This  behavior  is  often  encountered  when  the  function  /  is  obtained,  in 
whole  or  in  part,  by  table  interpolation.  With  a  better  initial  guess,  the 
method  would  have  succeeded. 


By  the  Newton-Raphson  formula. 

f(xi) 

x,+'~* ’  fVi)  ’ 

(04) 

50  that  /(*.) 

61+1  *  f'(xi)  • 

(05) 

When  a  trial  solution  X{  differs  from  the  true  root  by  e;,  we  can  use  (G3)  to  express  /(*{), 
f'(xi)  in  (G4)  in  terms  of  e2  and  derivates  at  the  root  itself.  The  result  is  a  recurrence 
relation  for  the  deviations  of  the  trial  solutions, 


€f+l  =  ~€i 


/"(*) 


(G6) 


2 /'(*) 

Equation  (G6)  says  that  Newton-Raphson  converges  quadratically.  Near  a  root,  the  number 
of  significant  digits  approximately  doubles  with  each  step.  This  very  strong  convergence 
property  makes  Newton-Raphson  the  method  of  choice  for  any  function  whose  derivative 
can  be  evaluated  efficiently  and  is  continuous  in  the  neighborhood  of  a  root,  and  we  have 
found  it  very  effective  in  correcting  locus  errors  A*  (see  Section  3.4). 


Even  when  Newton-Raphson  is  rejected  for  the  early  stages  of  convergence  (because 
of  its  poor  global  convergence  properties),  it  is  very  common  to  “polish  up”  a  root  with 
one  or  two  steps  of  Newton-Raphson,  which  can  multiply  by  two  or  four  its  number  of 
significant  figures! 


For  an  efficient  realization  of  Newton-Raphson  the  user  provides  a  function  subroutine 
which  evaluates  both  f(x)  and  its  first  derivative  f'(x)  at  the  point  x.  The  Newton- 
Raphson  formula  can  also  be  applied  using  a  numerical  difference  to  approximate  the  true 
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local  derivative. 


_  (G7) 
dx 

This  is  not,  however,  a  recommended  procedure  for  the  following  reasons:  (i)  You  are 
doing  two  function  evaluations  per  step,  so  at  best  the  superlinear  order  of  convergence 
will  be  only  \/2.  (u)  If  you  take  dx  too  small  you  will  be  wiped  out  by  roundoff,  while 
if  you  take  it  too  large  your  order  of  convergence  will  be  only  linear,  no  better  than  using 
the  initial  evaluation  f\x o)  for  all  subsequent  steps.  Therefore,  Newton-Raphson  with 
numerical  derivates  is  (in  one  dimension)  always  dominated  by  the  secant  method.  (In 
multidimensions,  where  there  is  a  paucity  of  available  methods,  Newton-Raphson  with 
numerical  derivates  must  be  taken  more  seriously.) 

The  following  subroutine  calls  a  user-supplied  subroutine  FUNCD  (X,  FN,  DF)  which 
returns  the  function  value  as  FN  and  the  derivative  as  DF.  We  have  included  input  bounds 
on  the  root  simply  to  be  consistent  with  previous  root-finding  routines:  Newton  does  not 
adjust  bounds,  and  works  only  on  local  information  at  the  point  X.  The  bounds  are  only 
used  to  pick  the  midpoint  as  the  first  guess,  and  to  reject  the  solution  if  it  wanders  outside 
of  the  bounds. 

FUNCTION  RTNEWT (FUNCD , XI , X2 , XACC) 

Using  a  the  Newton-Raphson  method,  find  the  root  of  a  function  known  to  lie  in 
the  interval  X1-X2.  The  root  RTNEWT  will  be  refined  until  its  accuracy  is  known 
within  ±XACC.  FUNCD  is  a  user- supplied  subroutine  that  returns  both  the  function 
value  and  the  first  derivative  of  the  function  at  the  point  X. 


P  ARAMS  TER  ( JMAX= 20)  Set  to  maximum  number  of  iterations . 

RTNEWT*  .5*  (X1+X2 )  Initial  guess. 

DOim  J=l,  JMAX 

CALL  7UNCD (RTNEWT, 7, DT) 

DX=F /DF 

RTNEWT =RTNEWT-DX 

IF ( (XI -RTNEWT) * (RTNEWT-X2) . LT . 0 . ) PAUSE  ’jumped  out  of  bracket# ' 
IF  (ABS  (DX)  .LT.XACC)  RETURN  Convergence. 
flU  CONTINUE 

PAUSE  ’RTNEWT  exceeding  maximum  iterations' 

END 


While  Newton-Raphson’s  global  convergence  properties  are  poor,  it  is  fairly  easy  to 
design  a  fail-safe  routine  that  utilizes  a  combination  of  bisection  and  Newton-Raphson.  The 
hybrid  algorithm  takes  a  bisection  step  whenever  Newton-Raphson  would  take  the  solution 
out  of  bounds,  or  whenever  Newton-Raphson  is  not  reducing  the  size  of  the  brackets  rapidly 
enough. 
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FUNCTION  RTSAFE(FUNCD , XI , X2 , XACC) 


Using  a  combination  of  Newton-Raphson  and  bisection,  find  the  root  of  a  function 
bracketed  between  XI  and  X2.  The  root,  returned  as  the  function  value  RTSAFE, 
will  be  refined  until  its  accuracy  is  known  within  ±XACC.  FUNCD  is  a  user  supplied 
subroutine  which  returns  both  the  function  value  and  the  first  derivative  of  the 
function. 


PARAMETER  ( MAXI  T= 1 0  0 )  Maximum  allowed  number  of  iterations. 

CALL  FUNCD (XI , FL , DF) 

CALL  FUNCD (X2 , FH , DF) 

IF (FL*FH .GE . 0 . )  PAUSE  ‘root  must  be  bracketed* 

IF  (FL .  LT .  0  . )  THEN  Orient  the  search  so  that  f(XL)  <  0. 

XL  =  X1 
XH=X2 


ELSE 

XH  =  X1 


XL  =  X2 
SWAP=FL 
FL  =  FH 
FH=SWAP 
END  I F 

RTSAFE= .5* (X1+X2) 
DXOLD=ABS (X2-X1) 

DX  =  DX  0 LD 

CALL  FUNCD (RTSAFE, F,DF) 
DO  [O  J=1 , MAXIT 


Initialize  the  guess  for  root, 
the  "step-size  before  last," 
and  the  last  step. 

Loop  over  allowed  iterations 


IF  (  (  (RTSAFE-XH)  *DF-F)  *  (  (RTSAFE-XL)  *DF-F)  .  GE  .  0  .  Bisect  if  Newton  out  of  range, 
.  OR .  ABS  (2 .  *F)  .GT .  ABS  ( DXOLD  *DF)  )  THEN  or  not  decreasing  fast  enough. 

DXO  LD  =  DX 
DX=  0 .5* (XH-XL) 

RTSAFE*XL+DX 
IF ( XL . EQ . RTSAFE )  RETURN 

ELSE 

DXOLD=DX 
DX=F / DF 
TEMP = RTSAFE 
RTSAFE=RTSAFE-DX 
IF (TEMP .EQ. RTSAFE)  RETURN 
END  IF 

IF (ABS (DX) . LT . XACC)  RETURN 
CALL  FUNCD (RTSAFE , F , DF) 

IF (F . LT . 0 . )  THEN 


Change  in  root  is  negligible. 
Newton  step  acceptable.  Take  it 


Convergence  criterion. 

The  one  new  function  evaluation  per  iteration. 
Maintain  the  bracket  on  the  root 


XLsRTSAFE 

FL=F 


ELSE 


XH=RTS AFE 
FH=F 


NDIF 

CONTINUE 


PAUSE  'RTSAFE  exceeding 

RETURN 

END 


maximum 


iterations  * 


For  many  functions  the  derivative  f'(x)  often  converges  to  machine  accuracy  before 
the  function  /(a?)  itself  does.  When  that  is  the  case  one  need  not  subsequently  update 
f'(x).  This  shortcut  is  recommended  only  when  you  confidently  understand  the  generic 
behavior  of  your  function,  but  it  speeds  computations  when  the  derivative  calculation  is 
laborious.  (Formally  this  makes  the  convergence  only  linear,  but  if  the  derivative  isn’t 
changing  anyway,  you  can  do  no  better.)  For  further  reading  on  this  subject,  see  [G1-G3]. 
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G.2.  Minimizing  a  Function  f(x) 

Here  we  employ  two  recipes,  DBRENT  [G3],  and  BRENT  [G3].  The  first  is  used  to 
compute  the  minimum  min^M{|Ar/(<^]i;f  )|},  whose  derivatives  were  analytically  derived  in 
Appendix  E.  The  second  is  used  to  compute  the  value  twait  of  fwa;t  for  which  the  miss 
distance  MD[tyja it)  reaches  a  minimum.  Here  we  do  not  have  analytically  determinable 
derivatives,  and  BRENT  computes  these  as  well  as  t^ait 


G.3.  Parabolic  Interpolation  and  Brent’s  Method  in  One-Dimension 

A  golden  section  search  is  designed  to  handle,  in  effect,  the  worst,  possible  case  of 
function  minimization,  with  the  uncooperative  linnimum  hunted  down  and  cornered  like  a 
scared  rabbit.  But  why  assume  the  worst?  If  the  function  is  nicely  parabolic  near  to  the 
minimum — surely  the  generic  case  for  sufficiently  smooth  functions  —then  the  parabola 
fitted  through  any  three  points  ought  to  take  Us  in  a  single  leap  to  the  minimum,  or  at 
least  very  near  to  it  (see  Figure  G.4).  Since  we  want  to  find  an  abscissa  rather  than  an 
ordinate,  the  procedure  is  technically  called  inverse  parabolic  interpolation. 


parabola  through  l,  2,  z 
parabola  through  l,  2,  4 


Figure  G.4.  Convergence  to  a  minimum  by  inverse  parabolic  interpolation.  A  parabola 
(dashed  line)  is  drawn  through  the  three  original  points  1,  2,  3  on  the 
given  function  (solid  line).  The  function  is  evaluated  at  the  parabola’s 
minimum,  4,  which  replaces  point  3.  A  new  parabola  (dotted  fine)  is 
drawn  through  points  1,  4,  2.  The  minimum  of  this  parabola  is  at  5, 
which  is  close  to  the  minimum  of  the  function. 
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The  formula  for  the  abscissa  x  which  is  the  minimum  of  a  parabola  through  three 
points  /(a),  /(&),  and  /(c)  is 

c)2[f(b)  ~  f(a)} 


,1  (6-a)2[/(fe)-/(c)]-(6 _ 

2  (b  -  a)[f(b)  -  /(c)]  -  (b  -  c)[f(b)  -  /(a)] 


(G8) 


as  you  can  easily  derive.  This  formula  fails  only  if  the  three  points  are  collinear,  in  which 
case  the  denominator  is  zero  (minimum  of  the  parabola  is  infinitely  far  away).  Note, 
however,  that  (G8)  is  as  happy  jumping  to  a  parabolic  maximum  as  to  a  minimum.  No 
minimization  scheme  that  depends  solely  on  (G8)  is  likely  to  succeed  in  practice. 


The  exacting  task  is  to  invent  a  scheme  which  relies  on  a  sure-but-slow  technique, 
like  golden  section  search,  when  the  function  is  not  cooperative,  but  which  switches  over 
to  (G8)  when  the  function  allows.  The  task  is  nontrivial  for  several  reasons,  including 
these:  ( i )  The  housekeeping  needed  to  avoid  unnecessary  function  evaluations  in  switching 
between  the  two  methods  can  be  complicated,  (it)  Careful  attention  must  be  given  to 
the  “endgame,”  where  the  function  is  being  evaluated  very  near  to  the  roundoff  limit  of 
equation  (G8).  (in)  The  scheme  for  detecting  a  cooperative  versus  noncooperative  function 
must  be  robust. 

Brent’s  method  [G4]  is  up  to  the  task  in  all  particulars.  At  any  particular  stage,  it 
is  keeping  track  of  six  function  points  (not  necessarily  all  distinct),  o,  b,  u,  v,  w,  and  x, 
defined  as  follows:  the  minimum  is  bracketed  between  a  and  b;  x  is  the  point  with  the 
very  least  function  value  found  so  far  (or  the  most  recent  one  in  case  of  a  tie);  w  is  the 
point  with  the  second  least  function  value;  v  is  the  previous  value  of  w;  u  is  the  point  at 
which  the  function  was  evaluated  most  recently.  Also  appearing  in  the  algorithm  is  the 
point  xm,  the  midpoint  between  a  and  b;  however,  the  function  is  not  evaluated  there. 

You  can  read  the  code  below  to  understand  the  method’s  logical  organization.  Men¬ 
tion  of  a  few  general  principles  here  may,  however,  be  helpful:  Parabolic  interpolation  is 
attempted,  fitting  through  the  points  x,  v,  and  w.  To  be  acceptable,  the  parabolic  step 
must  (i)  fall  within  the  bounding  interval  (a,  b),  and  (it)  imply  a  movement  from  the  best 
current  value  x  that  is  less  than  half  the  movement  of  the  step  before  last.  This  second  cri¬ 
terion  insures  that  the  parabolic  steps  are  actually  converging  to  something,  rather  than, 
say,  bouncing  around  in  some  nonconvergent  limit  cycle.  In  the  worst  possible  case,  where 
the  parabolic  steps  are  acceptable  but  useless,  the  method  will  approximately  alternate  be¬ 
tween  parabolic  steps  and  golden  sections,  converging  in  due  course  by  virtue  of  the  latter. 
The  reason  for  comparing  to  the  step  before  last  seems  essentially  heuristic:  experience 
shows  that  it  is  better  not  to  “punish”  the  algorithm  for  a  single  bad  step  if  it  can  make 
it  up  on  the  next  one. 


Another  principle  exemplified  in  the  code  is  never  to  evaluate  the  function  less  than 
a  distance  TOL  from  a  point,  already  evaluated  (or  from  a  known  bracketing  point).  The 
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reason  is  that,  as  we  saw  in  equation  (G8),  there  is  simply  no  information  content  in  doing 
so:  the  function  will  differ  from  the  value  already  evaluated  only  by  an  amount  of  order  the 
roundoff  error.  Therefore  in  the  code  below  you  will  find  several  tests  and  modifications 
of  a  potential  new  point,  imposing  this  restriction.  This  restriction  also  interacts  subtly 
with  the  test  for  “doneness,”  which  the  method  takes  into  account. 

A  typical  ending  configuration  for  Brent’s  method  is  that  a  and  b  are  2  x  x  x  TOL  apart, 
with  x  (the  best  abscissa)  at  the  midpoint  of  a  and  6,  and  therefore  fractionally  accurate 
to  ±T0L. 

Indulge  us  a  final  reminder  that  TOL  should  generally  be  no  smaller  than  the  square 
root  of  your  machine’s  floating  point  precision. 

FUNCTION  BRENT (AX , BX , CX , F , TOL , XMIN) 

Given  a  function  F,  and  given  a  bracketing  triplet  of  abscissas  AX,  BX,  CX  (such 
that  BX  is  between  AX  and  CX,  and  F(AX)  and  F(CX)),  this  routine  isolates  the 
minimum  to  a  fractional  precision  of  about  TOL  using  Brent’s  method.  The  abscissa 
of  the  minimum  is  returned  as  XMIN,  and  the  minimum  function  value  is  returned 
as  BRENT,  the  returned  function  value. 

For  further  reading  on  this  subject,  see  [G5]. 


G.4.  One-Dimensional  Search  with  First  Derivatives 

Here  we  want  to  accomplish  precisely  the  same  goal  as  in  the  previous  section,  namely 
to  isolate  a  functional  minimum  that  is  bracketed  by  the  triplet  of  abscissas  a  <  b  <  c, 
but  utilizing  an  additional  capability  to  compute  the  function’s  first,  derivative  as  well  as 
its  value. 

In  principle,  we  might  simply  search  for  a  zero  of  the  derivative,  ignoring  the  function 
value  information,  using  a  root  finder  like  RTFLSP  or  ZBRENT  [G3].  It  doesn’t  take  long  to 
reject  that  idea:  How  do  we  distinguish  maxima  from  minima?  Where  do  we  go  from  initial 
conditions  where  the  derivatives  on  one  or  both  of  the  outer  bracketing  points  indicate  that 
“downhill”  is  in  the  direction  out  of  the  bracketed  interval? 

We  don’t  want  to  give  up  our  strategy  of  maintaining  a  rigorous  bracket  on  the  min¬ 
imum  at  all  times.  The  only  way  to  keep  such  a  bracket  is  to  update  it  using  function 
(not  derivative)  information,  with  the  central  point  in  the  bracketing  triplet  always  that 
with  the  lowest  function  value.  Therefore  the  role  of  the  derivatives  can  only  be  to  help 
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PARAMETER  (ITMAX-100,  CGOLD- . 3819660,  ZEPS-1 . 0B-10) 

Maximum  allowed  number  of  iterations,  golden  ratio,  and  a  small  number  which  protects  against  trying  to  achieve  fractional 
accuracy  for  a  minimum  that  happens  to  be  exactly  zero. 


A-MIN(AX,CX)  A  and  B  must  be  in 

B-MAX(AX,CX)  though  the  input  ab 

^ = Initializations 

X«V 

FX-F(X)  This  will  be  the  dis 

FV-FX 

FW-FX 

D0[ll|  ITER-1, ITMAX  Main  program  loop. 

XM—  0.5* (A+B) 

T0L1— T0L*ABS (X) +ZEPS 
TOL2-2. *TOLl 

IF  (ABS  (X-XM)  .  LE  .  (T0L2- .  5*  (B-A)  )  )  GOTO  3  Test  for  done  here. 
IF<ABS(E)  .GT.TOL1)  THEN  Construct  a  trial  pai 

R-(X-M) * (FX-FV) 

Q-(X-V) * (FX-FW) 

P* (X-V) *Q— (X-W) *R 
0*2 . * (0-R) 

IF  (0  •  GT .  0  . )  P— P 
Q-ABS (Q) 

ETEMP-E 
X  -D 

IF (ABS (P) . GE .ABS ( . 5*Q*ETEMP)  .OR.  P.LE.Q*(A-X)  . 
P . GE . Q* (B— X) )  GOTO  1 

The  above  conditions  determine  the  acceptability  of  the  parabolic  fit  Here  it  is  o  Jc. 


A  and  B  must  be  in  ascending  order, 
though  the  input  abscissas  need  not  be. 
Initializations 


This  will  be  the  distance  moved  on  the  step  before  last 


Construct  a  trial  parabolic  fit 


P . LE  . 0*  (A-X)  .OR. 


U  *  X  +  D  Take  the  parabolic  step. 

IF (U-A.LT.TOL2  .OR.  B-U . LT.T0L2)  D-SIGN (TOL1 , XM-X) 


GOTO  2 
END  I  F 

1  IF(X.GE.XM)  THEN 

E-A-X 

ELSE 

E-B-X 
END  I  F 
D  — CGOLD  *E 

2  IF (ABS (D) .GE.TOL1)  THEN 

U-X+D 

ELSE 

U— X+SIGN (TOL1 , D ) 

END  IF 
FU—F (U) 

IF(FU.LE.FX)  THEN 

IF(U.GE.X)  THEN 
A-X 

ELSE 

B-X 

END  I  F 
V- W 
FV-FW 

w-x 

FW-FX 

X-U 

FX=FU 

ELSE 

IF(O.LT.X)  THEN 
A-U 

ELSE 

B-U 
END  I F 

IF (FU . LE . FW  .OR.  H.EQ.X)  THEN 
V  =  W 
FV-FH 
W-U 
FW=FU 

ELSE  IF  (FU .  XiE  .FV  .OR.  V.EQ.X  .OR. 
V-U 
FV-FU 
END  IF 
END  IF 

flH  CONTINUE 

PAUSE  'Br«nt  •xc««d  maximum  itarationa. > 

3  XM IN— X 
BRENT-FX 
RETURN 
END 


Skip  over  the  golden  section  step. 

We  arrive  here  for  a  golden  section  step, 
which  we  take  into  the  larger  of  the  two  segments. 


Take  the  golden  section  step. 

Arrive  here  with  D  computed  either  from  parabolic  fit, 
or  else  from  golden  section. 


This  is  the  one  function  evaluation  per  iteration,  and  now 
we  have  to  decide  what  to  do  with  our  function  evaluation. 
Housekeeping  follows: 


V.EQ.W)  THEN 


Done  with  housekeeping.  Back  for  another  iteration. 


Arrive  here  ready  to  exit  with  best  values. 
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us  choose  new  trial  points  within  the  bracket. 

One  school  of  thought  is  to  “use  everything  you’ve  got”:  Compute  a  polynomial  of 
relatively  high  order  (cubic  or  above)  which  agrees  with  some  number  of  previous  function 
and  derivative  evaluations.  For  example,  there  is  a  unique  cubic  that  agrees  with  function 
and  derivative  at  two  points,  and  one  can  jump  to  the  interpolated  minimum  of  that  cubic 
(if  there  is  a  minimum  within  the  bracket).  Suggested  by  Davidon  and  others,  formulas 
for  this  tactic  are  given  in  Acton  [Gl] . 


We  like  to  be  more  conservative  than  this.  Once  superlinear  convergence  sets  in,  it 
hardly  matters  whether  its  order  is  moderately  lower  or  higher.  In  practical  problems  that 
we  have  met,  most  function  evaluations  are  spent  in  getting  globally  close  enough  to  the 
minimum  for  superlinear  convergence  to  commence.  So  we  are  more  worried  about  all  the 
funny  “stiff”  things  that  high  order  polynomials  can  do  and  about  their  sensitivities  to 
roundoff  error. 

This  leads  us  to  use  derivative  information  only  as  follows:  The  sign  of  the  derivative 
at  the  central  point  of  the  bracketing  triplet  a  <  b  <  c  indicates  uniquely  whether  the 
next  test  point  should  be  taken  in  the  interval  (a,  b)  or  in  the  interval  (fc,  c).  The  value 
of  this  derivative  and  of  the  derivative  at  the  second-best-so-far  point  are  extrapolated  to 
zero  by  the  secant  method  (inverse  linear  interpolation),  which  by  itself  is  superlinear  of 
order  1.618.  (The  golden  mean  again:  See  Acton,  p.  57  [Gl].)  We  impose  the  same  sort 
of  restrictions  on  this  new  trial  point  as  in  Brent’s  method.  If  the  trial  point  must  be 
rejected,  we  bisect  the  interval  under  scrutiny. 

Observe  that  we  have  emphasized  the  use  of  derivative  information  in  one-dimensional 
minimization.  But  beware  of  functions  whose  computed  “derivatives”  don’t  integrate  up 
to  the  function  value  and  don’t  accurately  point  the  way  to  the  minimum,  usually  because 
of  roundoff  errors,  sometimes  because  of  truncation  error  in  the  method  of  derivative 
evaluation. 

You  will  see  that  the  following  routine  is  closely  modeled  on  BRENT  in  the  previous 
section. 
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FUNCTION  DBRENT(AX,  BX,  CX,  F,  DF,  TOL,  XMIN) 

Given  a  function  F  and  its  derivative  function  DF,  and  given  a  bracketing  triplet  of 
abscissas  AX,  BX,  CX  (such  that  BX  is  between  AX  and  CX,  and  F(BX)  is  less  than 
both  F(AX)  and  F(CX)),  this  routine  isolates  the  minimum  to  a  fractional  precision 
of  about  TOL  using  a  modification  of  Brent’s  method  that  uses  derivatives.  The 
abscissa  of  the  minimum  is  returned  as  XMIN,  and  the  minimum  function  value  is 
returned  as  DBRENT,  the  returned  function  value. 


PARAMETER  (ITMAX-100,  ZEPS-1 . 0E-10) 

Comments  following  will  point  out  only  differences  from  the  routine  BRENT.  Read  that  routine  first 


LOGICAL  OKI,  OK2 
A-MIN (AX , CX) 
B-MAX(AX,CX) 

V-BX 

w-v 

x-v 

E-0 

fx-f (x) 

rv-rx 

fw-fx 

DX-Df (X) 

DV«DX 
DW— DX 

DO  EH  ITER-1, ITMAX 
XM— 0 . 5* (A+B) 


Will  be  used  as  flags  for  whether  proposed 
steps  are  acceptable  or  not 


All  our  housekeeping  chores  are  doubled  by  die  necessity  of 
moving  derivative  values  around  as  well  as  function  values. 


Initial i 7i»  these  D’s  to  an  out-of-bracket  value 


TOLl— TOL* ABS (X) +EEPS 
TOL2-2 . *TOLl 

IP(ABS(X-XM)  .LI.  (TOL2-  .  5*  (B-A)  )  )  GOTO  3 
IF (ABS (E) . GT . TOLl)  THEN 
Dl-2.* (B-A) 

gjSrS.'Sg 

Which  of  these  two  estimates  of  D  shall  we  take?  We  will  insist  that  they  be  within  die  bracket, 
and  on  the  side  pointed  to  by  the  derivative  at  X. 


Ul-X+Dl 

U2-X+D2 

OKI— ( (A-Ul) * (Ul-B) . GT . 0 . )  .AND. 
OK2— ( (A-U2) * (U2-B) . GT . 0 . )  .AND. 
OLDE-E 
E-D 

IF ( . NOT .  (OKI . OR . OK2  )  )  THEN 
GOTO  1 

ELSE  IF  (OKI  .AND.  OK2)  THEN 

IF (ABS (Dl) . LT .ABS (D2) )  THEN 


(DX*Dl . LE . 0 . ) 

(DX*D2 . LE . 0 . ) 

Movement  on  the  step  before  last 
Take  only  an  acceptable  D,  and  if  both  are 
acceptable,  then  take  the  smallest  one. 


D— Dl 

ELSE 


1 


2 


D-D  2 
ENDIF 

ELSE  IF  (OKI)  THEN 
D— Dl 

ELSE 

D-D  2 
ENDIF 

IF (ABS (D) .GT.ABS (0.5*OLDE) )  GOTO  1 
U-X  +  D 

IF (U-A.LT.TOL2  .OR.  B-U.LT.TOL2)  D-SIGN (TOLl , XM-X) 

GOTO  2 

IF  (Dx!ge  .  0 . )  THEN  Decide  which  segment  by  the  sign  of  the  derivative. 

E-A-X 

ELSE 

E-B-X 

n  -  n  1  *  Bisect,  not  golden  section. 

;;°AB3(D).<».T0L1)  THIN 
U-X  +  D 
FU-F (U) 
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ELSE 

U-X+SIGN (TOL1, D) 

ru-r (u) 

ir (ru.GZ.rx)  goto  3 
ENDir 
DU-DE (0) 

ir(ru.LE.rx)  them 

ir(U.GE.X)  THEM 
A-X 

ELSE 

B-X 
EMDir 

V -M 

rv-rw 

DV-DW 

W -X 

rw-rx 

DW-DX 
X-U 

rx-ru 

DX-DU 

ELSE 

ir(U.LT.X)  THEM 

A-U 

ELSE 

B-U 
END  I F 

ir(ru.LE.rw  .or.  w.eq.x)  them 
v- w 
rv-rw 

DV-DW 

w-u 

rw-ru 

DW-DU 

ELSE  ir(r0.LE.rV  .OR.  v.eq.x  .or.  V.EQ.W)  THEM 
v-u 

rv-ru 

DV-DU 

EMDir 

EMDir 

EIJ  CONTINUE 

PAUSE  1 DBRENT  •zc««d«d  maximum  it«rationa. 1 
3  XMIM— X 
DBREMT-rX 
RETURN 
END 


REFERENCES  AND  FURTHER  READING: 


[Gl]  Acton,  Forman  S.  1970,  Numerical  Methods  That  Work  (New  York:  Harper  and  Row), 
Chapter  2. 

[G2]  Ralston,  Anthony,  and  Rabinowitz,  Philip.  1978,  A  First  Course  in  Numerical  Analy¬ 
sis ,  2nd  ed.  (New  York:  McGraw-Hill). 

[G 3]  Ortega,  J.,  and  Rheinboldt,  W.  1970,  Iterative  Solution  of  Nonlinear  Equations  in 
Several  Variables  (New  York:  Academic  Press). 

[G4]  Brent,  Richard  P.  1973,  Algorithms  for  Minimization  without  Derivatives  (Englewood 
Cliffs,  N.  J.:  Prentice-Hall),  Chapter  5. 

[G5]  Forsythe,  George  E.,  Malcolm,  Michael  A.,  and  Moler,  Cleve  B.  1977,  Computer  Meth¬ 
ods  for  Mathematical  Computations  (Englewood  Cliffs,  N.J.:  Prentice-Hall  ),  §8.2. 


If  the  minimum  step  in  the  downhill  direction 
takes  us  uphill,  then  we  are  done. 

Now  all  the  housekeeping,  sigh. 
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APPENDIX  H 


Solving  for  tlj2 


We  solve 


Define  <fMC 


\JmccV°\  =  /  \Vm{t)\<It  for  tw 

M  Jo 

I Xw  —  Xc/2\.  There  are  three  cases  was  found  in  Chapter  3). 


for  ^u/2 


H.l.  Case  1.  X^f*ey  >  Hq ,  X*  ’ey  >  Ho • 


This  case  is  simple,  since  the  MaRV  speed  |VAf(T)|  —  I^mI>  a  constant. 


Hence 


,,/2  _  (Vc-x°„)-yA  _  ,ff/2 

-  /  I  t  rfl  l\0  LIC  ’ 


(IW 


/ .  ^t/2  x/2 

where  f »/?  is  the  solution  to  the  equation  |XC  —  =  f0IC  |V/(r)|dr,  and  Xc  was 

iL/  ir/2 

found  in  Chapter  3.  Note  that  tj  may  be  negative. 


H.2.  Case  2.  Xlf,ey  >  Ho ,  X^2  ’ey  <  Ho* 


The  MaRV  speed  |Vm(t)|  remains  constant  at  |V^r|  until  r  =  to  =  Hq  —  X°M'ey /V^ 'ty 
(Equation  (2.4)).  In  r0  seconds,  M  travels  |V^|r0  distance.  Hence,  all  we  need  is  to  solve 
the  following  equation  for  7: 

R=\X°m-  X;,2\  -  IVj&Iro  =  f  \VM(r)\ir  ,  (H2) 

Jo 

where 

|V„(T)|  -  llfil  [l  -  5^ 


We  obtain 


whose  solution  is: 


,7r/2  ,7t/2 

7  t JQ  *w  To 


rO  I  AyO  12 
mI-VI^mI  kDM 


\VM\/kDM 


—  kj DM  1  —  4  / 1  — 


IV^IfepM 
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or 


iv/2  = 


2  R 


\V°M\kDM 


+  To  “  / 


7r/2 

jc 


tt/2 

lw 


kDM 


1 


(H5) 


(H6) 


H.3.  Case  3.  <  H0l  X?/2-ey  <  H0. 

This  case  is  similar  to  H.2:  solve  the  following  equation  for  tw^  ■ 


The  solution  is 


l 


?r/2 

tic  +t 


x/2 


|Flf(r)|<fr 


(HT) 


.x/2 

f'W 


=  kf)]\f 


x/2 

l7C 


(H8) 
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Have 


Have  rj  (t. 


1.1. 


1.2. 


A  = 


B  = 


APPENDIX  I 

Some  Derivatives  for  the  Predictor-Corrector  Step  A* 


4>m)  —  ^(*)e2/M  5 

(11) 

|r^|  =  (^2  +  J52)1/2  . 

(12) 

,  At/,  4>*J  (<t>*M))  for  any 

tWi  4*j 

*7 

-i 

—A  sin  0®  +  B  cos  0® 

(13) 

=  cos 

c  1  c 

.  Ir/(^W7  At/,  <^f)|  _ 

A  (t, 

'!>■>  At/,  <f>M) 

=  k/HI 

-M*)l  7 

(14) 

A*(t, 

w,  At], 

=  min  {A(tw,  At/,  <£m)}  • 

<t>M 

(15) 

dA*  _  9\r'j\  _  d\rj I 
dAtj  dAt j  dAtj 


(16) 


Ml 

dAtj 


(A2  +  H2ri/2 


/  .  dA  D  dB  \ 

v Arni+Bm, )  ’ 


(17) 


m{4>*m)\  (sm<f>*M)  -  dMc(tw,  Atj)  -  dci(tw,  At/)  (cos  0°) 

-  |r/(Af/)|(sin0®)  ,  (18) 

m(<^m)I(cos 4>*m)  -  |rjif(Atjif)|  -  dci(tw ,  At/)  ( sin  0® )  +  |r/(At/)|  (cos0®)  ,  (19) 


l*if  -  *«<<.)!  - 

dMc{iw,  At/)  = 

/  |Viif(r)|dr  , 

Jo 

(110) 

1 

* 

1 

os 

»Ati 

dci{ty)i  At/)  = 

/  1  !"/(  t  )  |  At  , 

/o 

(Ill) 
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^KT,  ~  (tin#Jf)8|r"U")l  +  lF«(A<i  + '»)! 

+  (cos^)|V,(A<;)|  -  (sin»c°)a!-^Jl!  , 


(112) 


d^M^Af)!  lrM(*w,  4>*m)\  fdXM{tw,  Atj)-ey  dA hM u  A4  x  /T1oX 

“m7~  = - 4 - { - eM, - +  9a 77(  M)  ,(I13) 


^t-U»  H-  Atf 

Xm(1wi  Atj)  =  +  /  |Vm(t)|  dr  eyo  , 

Jo  M 

At/)  . 

- At/)|evo  , 


A<j,  m ) 
dA  tj 


' Rzmvmxm  (Ti  /^»  (sill  4>m) 

n  a  \  1  drM(tw,  < 

( 1  -  cos  0M)  eyMJ - 


drMjtw Atj) 
dAtj 


d\ri(Atj)\  dXj(Atj)-ey 

aiAij)  =  |r/(A(')l  -  9A(;ft0 


XI(At,)  =  X$  +  /  |V/(T)|*-eiUr,  , 

Jo 


8-Y;(A ij) 
dAtj 


=  mAtI)\eo.x,  , 


(114) 


(115) 


(116) 


(117) 


(118) 


(119) 


+  (cos,o)a^M  , 

A<j)  _  Afj)ffX m(*wi  Atjyey 

dA  tj  ~  h0  dAt! 


(120) 


(121) 


gN  _  rrOwi  A/j,  <ftj)  dXjjAtj) '  dAhijAti,  <f>j) 
dAtj  h0  dAtj  €y  dAtj 


(122) 
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AM  A </) 


h 

dAhj(Atj) 

dAtj 

dAhj(Atj) 
d\ri(Ati)  | 

g|r/(A<j)| 

0A</ 

dAhj(Atj) 

d<[>i 

dfa 

dAtj 

d<h 

d\r'j\ 


W(|r/|) 

«/(|rfr|) 

ML 

dAtj 

dB 


dfo 

dA 


^  Z  rur  ^ 

MAi/)|  (1  -  cos<^j)  - - - — +  (sin  0j)  eyo 

\ezm  X  eVf  I  J 


(123) 


=  COS 


-1 


5  cos  0®  —  A  sin  9 ® 


#A/ij(Aij)  9|rj(A//)|  ^  5A/*j(Atfj)  #</>/ 


#|r/(A<j)|  $A</  #A/j  ’ 


(124) 

(125) 


eZM  x  eFo 

(1  -  cos  <f>j )  - - 1-  (sin  <f>i)ev  o 

I  eZM  X  evy 


-J7  ) 


(126) 


was  developed  earlier  (Eq.  (D.19))  , 


(127) 


=  MAij)|e, 


€zm  X 


(sin0/)1 - —  +  (cos  <j>i)ev  o 

I  eZM  X  Cypl 


d<f>j  d\r'j\  d<fo  dB  d<j>j  dA 
d\r'j\  dAt 7  +  dAt I  +  dA  dAtj  ’ 


5 


=  ^7)  arccos(W(|^|))  =  - 


(l  —  77  (k/l)2) 1/2  ^ 


0  <  arccos  (U(\r\\))  <  ir  , 


B  cos  0®  —  A  sin  9® 


(128) 

(129) 


(130) 

(131) 


A  sin  0®  —  B  cos  0 ® 

TUJ\  2 


(132) 


was  derived  earlier  (Eqs.  (112),  (117),  (120))  ,  (133) 


cos  0 ® 

(l-M(|ri|)2)1/2|ri 

jg(g) 

(l  —  77  (lr/l)2) 1/2  lrjl 


(134) 


(135) 
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dB 

dAtj 


derived  earlier  (Eqs.  (120))  , 


(136) 


dA 

dAtj 


derived  earlier  (Eqs.  (112)) 


(137) 


Correction  Step: 


Choose  correction  step: 


Error: 


A*  = 


dA* 

dAtj 


At j 


(138) 

(139) 
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List  of  Symbols 


a:  Angle  of  attack 

«/l  (g/2):  Acceleration  of  interceptor  of  type  1  (2) 

ay  fcth  asset 

Cjy.  Drag  coefficient 

Cj:  Instantaneous  center  of  maximal  interceptor  turn 

Cy.  Lift  coefficient 

Cm-  Instantaneous  center  of  maximal  MaRV  turn 

Cq:  Costs  incurred  by  allocating  defense  resources  to  0q 
C(Ia):  Lifecycle  cost  of  one  interceptor  of  type  a 
8:  Interceptor  and  defense  system  delay  time 

A:  Geometric  distance  between  MaRV  and  interceptor  loci  (trajectories) 

AA tj:  Size  of  A tj  correction  step  used  to  bring  MaRV  and  interceptor  loci 
together 

A  hi:  Change  in  altitude  of  interceptor  I  during  banking  displacement 
angle  (f>j 

Ah]\f.  Change  in  altitude  experienced  by  MaRV  when  it  turns  through  an 
angular  displacement  <f>j \j 

A R-Ml{iw,  A t°j):  Total  distance  traveled  by  both  vehicles  for  waiting  time  tw  and 
interceptor  flight  time  At®  (At®  +  tw  =  A t°M) 

At/:  Flight  time  of  interceptor  I 

At®:  Interceptor  flight  time  at  which  d°MI  is  measured  (flight  time  to  reach 
turning  point) 

Atji/:  Flight  time  of  MaRV  since  last  state  update 

dCi(At°M ):  Distance  from  interception  point  Xc  to  interceptor  turning  position 
(at  turning  time  A t°M  =  At®  +  tw) 

dPii/rn.  Distance  from  current  MaRV  position  to  interception  (“collision”) 

M°  point  Vc 

dMc{&t°M):  Distance  from  interception  point  Xc  to  MaRV  turning  position 
(at  turning  time  A t°M) 
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<p  • 

aMl' 


ev?  \ev°) 

ey 

F 

u 

(7,  fi,  0) 
Ho: 
Hko: 
Ik'- 
Iku  • 
K : 


K 

k Dll  (&DJ2) 
kDM 
KRaq 
L-. 


if 

x 

M 

AT,- 


MD. 


aq 

m 


712  aql 

maqs: 


Geometric,  estimate  of  MaRV-interceptor  distance  at  banking  (turning) 
time  A t°M  =  A t°j  +  (for  initialization  purposes) 

Unit  vector  along  Vq  (V^) 

Unit  vector  along  y-axis  (altitude) 

Vehicle  centripetal  force 

Footprint  of  Oq  (a  random  function) 

Rotation  angles  of  (xjm ,  zm )  coordinates 

Effective  height  of  atmosphere 

Keep-out  altitude 

Interceptor  pool  associated  with  asset  a*, 
uth  interceptor  in  interceptor  pool  I 
Collection  of  assets 
Quantity  of  assets 

Deceleration  constant  of  interceptor  of  type  1  (2) 

MaRV  deceleration  constant 

Kill  radius  of  an  interceptor  of  type  a  against  threat  0q 
Lift  force 

Loss  incurred  if  object  0q  leaks  through 
Leakage  loss  incurred  if  Oq  is  a  BRV 
Leakage  loss  incurred  if  0q  is  a  MaRV 
Collection  of  MaRVs 
Mach  number 
jth  MaRV 

Miss  distance  of  interceptor  of  type  a  against  threat  0q 
Vehicle  mass 

Quantity  of  defense  resources  of  type  a  allocated  to  0q 
Number  of  interceptors  of  type  a  launched  from  site  s  towards  0q 
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m 


aq 

mT: 


m°n  (m/2) 


m 


M  ' 


Oq\ 

<t>r 

rr- 


<f>M 


*  . 
M  ' 


P 


Kq 


^Kq 

P  Lq 
PB 

yLq 

PM 

P  Lqs 


nD. 

P<1  ■ 
pf: 
Pf: 


Pqk  (Pfk): 
Q- 
Qi' 


Q*: 

P- 

Rq: 


Optimal  value  of  variable  maq 
Mass  of  interceptor  I 

Payload  (net  mass)  of  type  1  (type  2)  interceptor 
Mass  of  MaRV 
9th  threat  object 

Banking  (turning)  angle  (angular  displacement)  of  interceptor  I 

Interceptor  angular  displacement  from  start  of  turn  to  point  of  locus 
contact  ( Q *) 

Banking  or  turning  angle  (angular  displacement)  of  MaRV 

MaRV  angular  displacement  from  start  of  turn  to  point  of  locus 
contact  (£?*) 

Probability  that  Oq  will  be  killed  if  an  interceptor  is  launched  at  the 
optimal  time  t^ait 

Probability  that  0q  will  be  killed  if  an  interceptor  is  launched  now 
Probability  that  Oq  leaks  through  defense 
Leakage  probability  if  0q  is  a  BRV 
Leakage  probability  if  0q  is  a  MaRV 

Probability  that  Oq  will  leak  through  if  a  single  interceptor  is  launched 
from  site  s 

Probability  that  0q  is  a  decoy 
Probability  that  Oq  is  a  MaRV 
Probability  that  Oq  is  a  BRV 

Probability  that  0q  is  headed  for  a if  0q  is  a  MaRV  (BRV) 

Total  quantity  of  threat  objects 

Index  of  the  last  threat  object  on  the  priority  list  for  which  a  sufficient 
quantity  of  interceptors  is  available 

Point  of  locus  contact 

Atmospheric  density 

Extremal  region 
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Ri 

R2 

r3 

0 

nMI 


A,  <i»  /w  i4  • 
Tfof  *AT 


*m(AV): 


M<)b 

5: 
c  • 

sign: 

r: 


< 


T0F: 


t*: 


tF- 

l0  . 


tB: 

BM  (*£/): 


(<})= 


Exoat mospheric  region 
Intermediate  region 
Tailchase  region 

Estimate  of  total  distance  separating  vehicles  at  most  recent  update 
time  to 

MaRV  rotation  matrix  between  action  coordinates  and  MaRV  coordinates 

Vector  collinear  with  r*j  from  interceptor  turn  center  to  a  point  on 
the  interceptor  locus 

Vector  from  interceptor  turn  center  to  point  Q  on  MaRV  locus 

Minimum  curvature  radius  vector  of  MaRV  at  t  =  A Im 

Threat  ranking  ratio  for  object  0q 

Minimum  radius  of  curvature  of  interceptor  I 

Curvature  radius  of  MaRV  at  time  t 

Effective  vehicle  area 

Collection  of  sites  each  of  which  contains  at  least  one  interceptor  of 
type  a 

Signum  function 

Clock  increment  (clock  “step  time”) 

Time  remaining  before  reaching  Hq 
Time  remaining  for  Oq  to  reach  H^0 

Incidence  (collision)  angle  of  interceptor  relative  to  MaRV 

Total  flying  time  of  threat  object  before  ground  impact 

Time  at  which  interceptor  reaches  Xj 

Flying  time  from  Ho  to  ground  impact 

Time  when  both  vehicles  start  banking 

Turning  time  required  by  MaRV  (interceptor)  to  reach  Q* 

Flying  time  of  MaRV  from  point  below  Ho  to  ground  impact 
MaRV  (interceptor)  turning  time  to  reach  Q* 


134 


twi  ^wait* 

^wait: 

t°‘ 

lw 

Ua: 
V : 
V: 
V?: 
Vn  (Vn): 

v?(vg,)- 
I ',r(v&Y. 

in* 

Wi-. 

Y*/2. 

J\.r 


Xjc: 

xp 

xb 

Xn  (Xn): 
XM(t): 

x*M(n: 

yO. 

Ag- 

(XM1  Vmi  ZJlf): 

(XM>  zm): 


Waiting  time 

Optimal  launch  delay  (waiting)  time 
Initial  waiting  time  (usually  t ^  =  0) 

Inventory  size  of  type  a  resources 

Value  function  of  assets  (usage  clear  from  context) 

Vehicle  velocity  (usage  clear  from  context) 

Launch  velocity  of  interceptor 
Velocity  of  interceptor  of  type  1  (2) 

Current  velocity  of  Oq  (MaRV) 

Velocity  of  Oq  (MaRV)  at  impact, 
vehicle  speed 

Opportunity  window  i  (i  =  1,2,  3) 

Point  on  the  extended  MaRV  trajectory  closest  to  the  interceptor 
launch  location 

Interception  point  in  the  tailchase  region 
Position  of  interceptor  when  reaching  Q*  (Xj  =  Q*) 

Launch  location  of  interceptor  I 
Displacement  of  interceptor  of  type  1  (2) 

Position  of  MaRV  at  time  i 
Position  of  MaRV  when  Xj  =  Q* 

Current  position  of  Oq 
MaRV  coordinates  (M- system) 

Action  coordinate  system  (A-system) 
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